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ABSTRACT 


To  meet  demands  for  limited  water  resources,  managers  of  water  supply  systems 
are  considering  innovative  management  methods.  "Conjunctive  use"  the 
coordinated  management  of  surface  water  and  groundwater — is  frequently 
espoused  as  a  method  for  addressing  these  demands  by  improving  the  reliability 
and  efficiency  of  water  management  activities. 

Water  management  is  often  heuristic,  based  on  prior  experience  managing 
existing  systems.  Because  conjunctive-use  is  relatively  new,  we  do  not  have 
heuristic  rules  for  managing  the  allocation  of  water  between  surface  water 
reservoirs  and  aquifers. 

This  thesis  presents  an  analytic  procedure  to  derive  rules  for  managing  a 
conjunctive-use  system.  In  this  thesis,  system  management  refers  to  the 
specification  of  optimal  controls:  decisions  regarding  the  release  and  allocation 
of  water.  The  procedure  is  applied  to  a  proposed  conjunctive-use  system  of 
practical  interest,  where  the  facilities  have  pre-determined  characteristics  and 
where  control  decisions  are  affected  by  uncertain  autocorrelated  inputs.  Optimal 
decisions  are  determined  by  applying  the  optimization  method  of  gradient 
dynamic  programming  to  a  numerical  model  of  the  conjunctive-use  system. 

The  resulting  optimal  decisions  demonstrate  that  application  of  the  analytic 
procedure  can  greatly  improve  the  management  of  a  conjunctive-use  system. 
These  decisions  are  significantly  more  efficient  than  those  that  result  from 
application  of  heuristic  rules.  Also,  the  analytic  procedure  incorporates  the  effect 
of  autocorrelation,  and  the  resulting  optimal  decisions  demonstrate  that 
autocorrelation  is  important  in  the  control  of  conjunctive-use  systems  because  of 
the  different  capabilities  and  constraints  that  surface  reservoir  storage  and 
aquifer  storage  present. 
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I.  Introduction 


Increasing  demands  on  limited  water  resources  are  leading  managers  of 
water  supply  systems  to  search  for  system  designs  and  control  methods  that  are 
efficient  and  innovative.  "Conjunctive  use"— the  coordinated  management  of 
surface  water  and  groundwater— has  received  much  attention  as  an  affordable 
and  environmentally  sound  method  for  enhancing  water  supply  systems. 
Increasingly,  conjunctive  use  is  proposed  as  a  solution  to  water  supply  problems 
(Sax  et  al,  1991;  EBMUD,  1992;  Western  Water,  1993),  and  it  is  likely  that  we  soon 
will  see  implementation  of  large-scale  conjunctive-use  projects.  As  yet,  however, 
there  have  been  few  projects  implemented  and  few  guidelines  exist  for  their 
management. 

Because  of  the  lack  of  heuristic  or  analytic  guidelines,  system  managers 
cannot  anticipate  how  they  should  optimally  control  new  or  modified  systems 
that  include  both  surface  water  and  groundwater.  In  particular,  the  different 
capabilities  and  constraints  of  surface-water-reservoir  (hereafter  abbreviated  to 
"reservoir")  and  aquifer  storage  result  in  optimal  controls  that  cannot  generally 
be  found  by  common  sense  or  heuristic  approaches.  A  purpose  of  this  thesis  is  to 
illustrate  how  we  can  find  optimal  controls  and  how  the  different  capabilities 
and  constraints  of  reservoir  and  aquifer  storage  affect  optimal  control.  I  will 
accomplish  this  by  modeling  and  solving  the  control  decisions  that  are  optimal 
for  a  conjunctive-use  system. 

A.  MOTIVATION  FOR  THIS  WORK: 

PLANNED  CONJUNCTIVE  USE  BY  A  WATER  SUPPLY  AGENCY 

Managers  of  water  supply  systems  often  hold  two  extreme  though 
conflicting  views  when  contemplating  the  addition  of  aquifer  storage  to  an 
existing  reservoir  system.  On  one  hand,  they  often  desire  to  view  aquifer  storage 
in  terms  of  an  equivalent  reservoir.  On  the  other  hand,  as  Lettenmaier  and  Burges 
[1979]  note,  they  tend  to  view  aquifer  storage  as  a  back-up,  to  be  used  in  times  of 
plenty  only  after  reservoirs  are  full  and  in  times  of  shortage  only  after  reservoirs 
are  empty.  Both  extreme  views  fail  to  recognize  the  different  capabilities  and 
limitations  of  reservoir  and  aquifer  storage. 

In  reality,  optimal  control  of  conjunctive-use  systems  is  far  from  either  of 
these  extreme  views.  This  thesis  draws  attention  to  some  of  the  fundamental 
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differences  between  aquifer  and  reservoir  storage,  as  well  as  the  benefits  of 
integrated  management.  If  we  are  to  operate  a  system  efficiently  and  reliably, 
management  rules  cannot  be  based  on  notions  that  groundwater  provides  only  a 
backup  supply  or  that  aquifer  storage  is  equivalent  to  reservoir  storage. 

I  base  my  comparison  of  the  two  storage  mechanisms  on  the  analysis  of  a 
hypothetical  conjunctive-use  system.  This  hypothetical  system  is  derived  from 
the  East  Bay  Municipal  Utilities  District  (EBMUD)  system  in  California.  EBMUD 
supplies  water  to  many  communities  along  the  eastern  shore  of  the  central  and 
southern  San  Francisco  Bay.  EBMUD  anticipates  a  growing  demand  for  water 
resources  in  its  district;  and  it  has  been  studying  ways  of  improving  its  ability  to 
meet  these  demands. 

Until  recently,  large  water  supply  agencies  such  as  EBMUD  relied  on  the 
expansion  of  reservoir  storage  to  meet  their  growing  needs.  Additional  storage 
allows  agencies  to  make  more  effective  use  of  fluctuating  and  uncertain  stream 
flows.  To  make  effective  use  of  its  fluctuating  and  uncertain  water  source  of  the 
Mokelumne  River,  EBMUD  initially  planned  to  develop  the  Buckhom  Reservoir. 
This  reservoir  would  have  added  capacity  to  an  extensive  reservoir  system 
consisting  of  five  reservoirs  in  EBMUD’s  service  district  and  two  large  reservoirs 
in  the  Sierra  Nevada  Mountain  Range. 

However,  two  events  have  thwarted  EBMUD’s  efforts  to  develop 
Buckhom  Reservoir.  The  first  was  a  set  of  legal  suits  brought  against  EBMUD  in 
1989.  Five  groups  charged  the  District  with  failure  to  consider  adequately  other 
options:  federal  and  state  law  requires  consideration  of  impacts  and  alternatives 
in  an  Environmental  Impact  Statement  /Report  (EIS/EIR)  before  construction  of 
a  Reservoir.  The  second  was  the  1992  elections  for  EBMUD's  board  of  directors. 
The  elections  produced  a  public  debate  on  the  construction  of  Buckhorn 
Reservoir  and  resulted  in  the  replacement  of  four  of  the  five  board  members  with 
individuals  who  campaigned  against  the  Reservoir. 

As  a  result,  Buckhorn  Reservoir  is  no  longer  a  serious  option  for 
enhancing  the  EBMUD  system.  Conjunctive  use,  on  the  other  hand,  has  become 
a  serious  option;  and  the  revised  EIS/EIR  identifies  it  as  the  "environmentally 
superior  alternative"  among  the  methods  considered  [EBMUD,  1992].  In  1993, 
the  EBMUD  board  of  directors  endorsed  conjunctive  use  as  the  primary  method 
for  enhancing  its  system. 

EBMUD’s  study  of  conjunctive  use  in  its  revised  EIS/EIR  originally 
motivated  this  work.  Though  the  District  initially  expressed  little  desire  to 
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implement  a  conjunctive-use  program,  it  has  now  become  the  primary  method 
for  enhancing  the  future  reliability  and  yield  of  its  water  supply.  It  is  likely  that 
other  utilities  facing  needs  similar  to  EBMUD  will  implement  conjunctive-use 
methods. 

Though  development  of  new  or  modified  conjunctive-use  systems  is  likely 
to  increase,  there  are,  at  present,  no  established  guidelines  that  can  aid  planners 
in  even  the  most  preliminary  evaluation  of  benefits.  I  intend  this  work  to 
provide  a  procedure  for  efficient  control  of  a  conjunctive-use  system,  and  thus 
clarify  how  planners  and  managers  may  evaluate  and  operate  such  systems. 

B.  DEVELOPMENTS  IN  RESERVOIR  MANAGEMENT  AND 

CONJUNCTIVE  USE 

Much  of  our  understanding  and  practice  of  water  management  is  the 
result  of  many  years  experience  regulating  water  systems.  Based  on  this 
experience,  managers  have  developed  operating  rules  that  minimize  risk  and 
expected  costs  arising  from  water  shortages  and  floods,  and  that  maximize  the 
benefit  of  hydropower  generation. 

For  many  systems,  the  years  of  operating  experience  have  produced  rules 
that  are  highly  efficient;  thus  analytic  optimization  may  provide  little  additional 
efficiency.  For  example,  optimization  of  an  existing  hydropower  system  by 
Kelman  et  al.  [1990]  improved  efficiency  only  slightly  (though  slight 
improvements  may  amount  to  a  considerable  financial  benefit).  Pre-existing 
rules  for  this  system  already  were  highly  efficient  after  years  of  experience 
gained  by  trial-and-error.  However,  Kelman  et  al.  developed  their  rules  with 
much  less  effort  and  without  subjecting  the  system  to  sub-optimal  control  while 
gaining  operating  experience.  Moreover,  they  were  able  to  develop  a  variety  of 
other  rules  for  different  constraints— an  ability  that  is  valuable  for  system 
planning  when  constraints  are  flexible. 

Though  many  existing  systems  are  efficiently  managed,  in  recent  years 
two  developments  have  increased  the  need  for  research  in  quantitative  methods 
for  water  systems  management.  The  first  development  is  that  some  of  the 
objectives  of  reservoir  management  have  changed.  Many  system  managers  must 
now  include  the  "costs"  and  "benefits"— though  ill-defined— of  qualities  such  as 
impacts  on  fisheries,  riparian  ecosystems,  scenery,  recreation,  and  water  quality. 
Assuming  that  these  values  are  quantifiable,  system  managers  face  the  problem 
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of  incorporating  these  qualities  in  new  operating  rules.  The  second  development 
is  that  reservoir  construction  has  become  very  difficult.  There  are  fewer 
remaining  sites  for  reservoir  construction  and  the  public  is  less  accepting  of  this 
as  a  solution  for  addressing  society's  increasing  competition  for  limited  water 
supplies.  As  Lettenmaier  and  Burges  [1979]  point  out,  "The  best  available  sites  for 
surface  storage  facilities  have  been  used  and  environmental  considerations  have 
eliminated  some  remaining  potential  sites."  As  a  result,  new  and  innovative 
management  techniques  are  being  developed.  One  of  the  most  popular  is  the 
integration  of  surface  and  groundwater  resources  in  a  single  "conjunctive-use" 
management  system. 

The  concept  of  conjunctive  use  has  been  around  for  many  years.  The  1957 
California  Water  Plan  pointed  out  that  in  the  Central  Valley  underground  storage 
within  200  feet  of  land  surface  is  about  six  times  larger  than  feasible  surface 
storage.  Also,  research  applied  to  the  conjunctive-use  problem  began  years  ago: 
one  of  the  first  quantitative  descriptions  of  a  conjunctive-use  management 
problem  was  by  Buras  in  1963.  Although  the  potential  benefits  of  using  the  large 
capacities  of  aquifer  storage  have  long  been  recognized,  few  projects  have  been 
undertaken  and  existing  literature  on  the  management  of  conjunctive-use 
systems  is  limited. 

Management  of  conjunctive-use  systems  presents  an  optimal  control 
problem  that  can  be  divided  into  two  phases:  modeling/simulation  and 
optimization.  While  the  modeling  and  simulation  of  water  supply  systems  is 
common,  optimization  is  used  infrequently.  Unfortunately,  available 
optimization  methods  have  been  limited  in  their  ability  to  accurately  address 
realistic  problems,  and  no  single  method  is  generally  applicable.  Because  of  this 
limitation,  most  problems  thus  far  have  been  simplified  excessively.  Burges  and 
Maknoon  [1975]  note  that,  "A  feature  of  nearly  all  the  literature  is  the  assumption 
that  one  or  several  parameters  or  variables  dominate  the  problem  at  hand." 
Consequently,  most  literature  provides  little  practical  guidance  to  system 
managers  seeking  to  plan  and  control  real  systems.  Not  only  have  past  efforts 
ignored  important  physical  characteristics  of  the  conjunctive-use  problem,  details 
about  important  non-physical  characteristics  arising  from  economic  and  legal 
factors  have  been  absent  or  lacking,  especially  when  compared  with  the  physical 
and  hydrologic  detail  normally  included  [Burges  and  Maknoon,  1975]. 
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C.  OBJECTIVES  OF  THIS  WORK: 


This  thesis  will  focus  on  the  effect  that  differences  in  storage  capacity  and 
transfer  rates  have  on  optimal  control  of  reservoir  and  aquifer  storage.  The 
storage  capacity  of  reservoirs  is  generally  more  limited  than  that  of  aquifers,  but 
the  transfer  rate  of  water  in  and  out  can  be  very  large.  In  contrast,  the  storage 
capacity  of  aquifers  is  generally  much  larger,  but  the  transfer  rate  is  limited. 

As  a  result,  efficient  conjunctive-management  practices  are  different  from 
management  practices  for  either  surface  or  subsurface  water  supplies  alone. 
Reservoir  storage  has  greater  ability  to  respond  to  needs  during  short-term  but 
severe  droughts,  while  aquifer  storage  has  greater  capacity  to  meet  water  needs 
during  long-term  but  moderate  droughts.  Efficient  management  of  a  system 
containing  both  reservoir  and  aquifer  storage  requires  the  optimal  use  of  the 
different  capabilities  and  limitations  of  each. 

This  work  presents  two  distinct  developments:  (1)  optimal  control  of  a 
conjunctive-use  system  that  does  not  rely  on  heuristic  groundwater  allocation 
rules  and  that  incorporates  the  influence  of  uncertain  autocorrelated  inputs,  and 
(2)  an  optimization  method  that  is  capable  of  solving  this  control  problem  while 
allowing  sufficient  complexity  to  provide  a  solution  of  practical  value.  Optimal 
control — release  and  allocation  decisions  that  best  achieve  the  objectives  of  a 
system — enables  us  to  observe  general  characteristics  of  control  decisions  aimed 
at  allocating  water  between  reservoir  and  aquifer  storage  and  at  hedging  against 
uncertainty. 

To  solve  a  conjunctive-use  problem,  I  apply  a  modified  gradient  dynamic 
programming  (GDP)  method.  GDP  is  a  dynamic  programming  method 
developed  by  Foufoula-Georgiou  and  Kitanidis  [1988]  to  solve  dynamic  control 
problems  of  greater  complexity  than  previously  possible.  In  this  thesis,  I  extend 
the  GDP  method  to  dynamic  control  problems  that  contain  autocorrelated 
stochastic  inputs.  An  autocorrelated  stochastic  input,  such  as  a  time-series  of 
stream  flows,  has  an  expected  future  value  that  depends  on  its  prior  values.  Our 
knowledge  of  these  prior  values  can  significantly  affect  our  control  decisions. 
Taking  advantage  of  GDP's  ability  to  solve  more  complex  problems— and  its 
extended  ability  to  include  autocorrelation— we  can  solve  conjunctive-use  control 
problems  with  fewer  simplifying  assumptions. 
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The  solution  of  the  conjunctive-use  control  problem  presented  in  this 
work  illustrates  two  important  observations  about  the  optimal  control  of  such 
systems.  The  first  observation  is  that  appropriate  control  of  aquifer  storage 
should  not  rely  on  heuristically  or  arbitrarily  defined  operating  rules.  This  is  a 
weakness  of  prior  research  on  optimal  control  of  conjunctive-use  systems:  prior 
research  generally  has  assigned  aquifer  storage  a  back-up  role  that  fails  to 
recognize  the  different  capabilities  and  limitations  of  reservoir  and  aquifer 
storage.  By  constraining  management  decisions,  this  back-up  role  has  resulted  in 
control  solutions  that  are  less  efficient  than  those  that  could  have  been  obtained 
by  the  approach  presented  here.  The  second  observation  is  that  autocorrelation 
can  significantly  affect  the  optimal  control  of  conjunctive-use  systems  because  of 
the  different  characteristics  of  reservoir  and  aquifer  storage. 

Prior  studies  of  conjunctive  use  systems  are  deficient  in  several  respects: 
(1)  they  lack  models  that  are  realistic  enough  to  solve  practical  management 
problems,  and  (2)  they  lack  an  optimization  approach  that  is  capable  of  solving 
any  but  the  most  simplified  problems.  The  two  developments  of  this  thesis 
partially  address  these  two  deficiencies.  Chapter  II  considers  characteristics  of 
optimally — vice  heuristically — derived  system  control  and  the  effect  of 
autocorrelation  on  optimal  control.  Chapter  IV  presents  an  extended  GDP 
method  that  can  allow  us  to  solve  optimal  control  problems  that  are  more 
complex  than  those  solvable  by  other  methods.  Chapters  III  considers  additional 
aspects  of  system  modeling,  and  Chapter  V  presents  conclusions,  practical 
extensions  of  this  work,  and  a  glossary  of  terms. 
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II.  Optimal  Control  of  Stochastic  Autocorrelated  Stream  Flows  in  a 
Conjunctively  Managed  Surface  Water  and  Groundwater  System 

A.  INTRODUCTION 

This  chapter  illustrates  and  solves  a  conjunctive-use  control  problem.  The 
conjunctive-use  model  is  based  on  an  existing  water  supply  system  that  currently 
uses  reservoirs  for  system  storage,  but  which  soon  may  include  aquifer  storage 
facilities.  The  optimal  control  solution  is  developed  without  the  use  of  heuristic 
control  rules  and  incorporates  the  effects  of  uncertain  autocorrelated  inputs. 

Through  this  example,  we  can  observe  significant  differences  in  the 
characteristics  of  reservoir  and  aquifer  storage  and  the  effect  that  these 
differences  have  on  optimal  control.  This  clarifies  the  benefits  of  conjunctively 
managing  surface  and  subsurface  water  resources,  thus  allowing  system 
planners  to  better  anticipate  the  costs  and  benefits  of  proposed  conjunctive-use 
systems.  Optimal  management  of  conjunctive-use  systems  takes  advantage  of 
the  strengths  of  both  reservoir  and  aquifer  storage  mechanisms,  achieving 
greater  system  efficiency  and  reliability  than  alternate  systems  using  either 
mechanism  in  isolation. 

In  this  section,  I  discusses  the  unique  character  of  the  conjunctive-use 
problem  by  focusing  on  the  different  capabilities  and  constraints  of  reservoir  and 
aquifer  storage.  In  addition,  I  discuss  the  influence  that  autocorrelation  can  have 
on  optimal  control,  an  influence  that  can  amplify  the  differences  between 
reservoir  and  aquifer  storage.  Such  influences  have  generally  been  ignored  by 
prior  conjunctive-use  research  as  reflected  in  the  literature  review  in  the  last  part 
of  this  section. 

The  following  sections  of  this  chapter  develop  and  solve  a  conjunctive-use 
problem.  Section  B  formulates  the  generic  control  problem  that  Section  C  applies 
to  a  specific  system  model.  Section  D  discusses  the  solution  method.  Sections  E 
and  F  present  results  and  conclusions  drawn  from  the  solution  of  the  specific 
problem. 

1.  Differences  Between  Reservoir  and  Aquifer  Storage 

In  a  dynamic  system,  current  decisions  influence  future  decisions  by 
altering  future  constraints  on  those  decisions  and  their  expected  cost.  As  a  result, 
system  operators  make  control  decisions  that  balance  current  performance 
against  future  performance.  When  operators  sacrifice  current  performance  for 
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improvements  in  the  future,  they  make  decisions  that  "hedge".  For  example,  a 
rational  operator  will  accept  rationing  of  current  water  supplies  if  the  cost  is 
more  than  offset  by  savings  from  less  severe  rationing  in  the  future. 

Uncertainty  plays  a  crucial  role  in  producing  control  decisions  in  dynamic 
system.  When  control  decisions  are  regulated  by  uncertain  inputs  (such  as  from 
stream  flow  or  imperfect  measurements  of  the  state  of  a  system),  it  is  not  possible 
to  identify,  in  advance,  a  single  control  schedule  that  optimally  balances  current 
and  future  performance  through  time.  Instead,  system  operators  must  determine 
control  decisions  in  "real  time" — they  must  continuously  update  control  decisions 
based  on  realized  and  predicted  values  of  inputs.  Because  realized  and  predicted 
values  of  inputs  improve  estimates  of  future  system  performance,  operators  can 
improve  control  decisions  that  balance  current  performance  against  expected 
future  performance. 

Stream  flow  and  many  other  hydrologic  phenomena  are,  by  their  nature, 
uncertain.  Conjunctive-use  systems  will  generally  depend  on  a  stream-flow 
source;  and,  as  a  result,  the  uncertainty  inherent  in  this  source  will  result  in 
optimal  control  decisions  that  hedge  supply  and  allocation  decisions.  However, 
unlike  control  systems  consisting  only  of  surface  water  reservoirs,  the  mix  of 
storage  mechanisms  in  conjunctive-use  systems  results  in  less  obvious  hedging 
policies.  The  different  constraints  on  the  use  of  surface  and  subsurface  storage 
suggest  that  the  optimal  control  of  the  different  components  may  be  quite 
different. 

Two  fundamental  differences  between  subsurface  and  reservoir  storage 
are  their  potential  size  and  rate  of  response.  Available  subsurface  storage 
frequently  exceeds  available  surface  storage  by  several  orders  of  magnitude.  On 
the  other  hand,  filling  and  draining  of  reservoirs  can  occur  at  rapid  rates,  while 
recharge  and  pumping  of  groundwater  can  be  significantly  constrained.  These 
fundamental  differences  indicate  that  optimal  control  of  reservoirs  and  aquifers 
should  be  quite  different:  while  we  might  store  water  in  reservoirs  for  use 
during  the  next  season  or  year,  we  might  store  water  in  the  subsurface  for  a 
drought  several  years  in  the  future.  Lettenmaier  and  Burges  [1979]  compare  the 
performance  of  groundwater  and  surface  water  systems  by  stating,  "In  contrast 
to  the  rather  long-term  failure  modes  encountered  in  excessive  reliance  on 
groundwater  supplies,  shorter  scale  (e.g.,  annual  or  seasonal)  failures  may  result 
from  exclusive  use  of  surface  supplies.  The  difference  in  time  scales  results 
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because  typical  surface  storage  reservoir  volumes  are  much  smaller  compared  to 
abstractions  than  are  groundwater  supplies." 

Because  of  these  differences,  water  supply  systems  should  be  designed  to 
employ  the  strengths  of  both  mechanisms,  resulting  in  greater  reliability  and 
lower  cost  than  a  system  using  either  method  in  isolation.  As  Burges  and 
Maknoon  [1975]  point  out,  "Whenever  multiple  sources  of  water  with  different 
characteristics,  as  is  the  case  with  groundwater  and  surface  water  systems,  are 
available,  it  may  be  possible  to  develop  an  operating  strategy  which  exploits  the 
different  characteristics  of  the  sources."  VJillis  and  Yeh  [1987]  recognized  that,  "By 
controlling  the  total  water  resources  of  a  region,  conjunctive  use  planning  can 
increase  the  efficiency,  reliability,  and  cost-effectiveness  of  water  use,  particularly 
in  river  basins  with  spatial  or  temporal  imbalances  in  water  demands  and  natural 
supplies." 

In  addition,  because  many  of  the  largest  water  supply  systems  have  relied 
entirely  on  reservoir  storage,  there  may  be  many  easy  opportunities  to  utilize 
aquifer  storage.  Initial  efforts  to  employ  aquifer  storage  may  be  far  cheaper  than 
expansion  of  reservoir  storage.  Lettenmaier  and  Burges  [1979]  found  that,  for  their 
model,  developing  aquifer  storage  to  buffer  against  variations  in  stream  flow  was 
about  an  order  of  magnitude  cheaper  than  developing  surface  storage. 

In  spite  of  the  potential  benefits,  few  water  supply  systems  have 
implemented  conjunctive-use  programs,  and  few  guidelines  exist  for  their 
operation.  Because  of  limited  experience  and  the  complexity  of  these  systems, 
efficient  control  of  new  or  modified  systems  is  difficult  to  determine  before  actual 
operation.  However,  development  of  heuristic  operating  rules  may  result  in 
unnecessarily  large  costs  early  in  a  project’s  life. 

2.  The  Importance  of  Autocorrelation  in  Optimal  Control 

Because  of  the  role  that  uncertainty  plays  in  producing  efficient  control 
decisions,  accurate  representation  of  this  uncertainty  is  essential.  One  should 
incorporate  any  information  that  can  improve  estimates  of  uncertain  inputs.  At 
time  scales  used  in  real-time  system  control,  autocorrelation  of  stream-flow  time 
series  is  high;  and,  as  a  result,  measurements  of  prior  stream-flows  are  essential 
to  accurate  estimation  of  future  flows.  Other  information  may  also  contribute  to 
estimates  of  future  stream  flow. 

Incorporating  autocorrelation  is  particularly  important  for  our 
understanding  of  the  conjunctive-use  problem.  Because  the  response 
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characteristics  of  reservoir  and  aquifer  storage  are  different,  optimal  control  of  a 
conjunctive-use  system  varies  with  the  degree  of  input  autocorrelation.  We  can 
consider  a  conjunctive-use  system  to  be  analogous  to  an  electrical  circuit  that 
consists  of  components  whose  impedances  (i.e.,  resistances  to  flows  of  electrons) 
depend  on  the  frequency  of  an  applied  signal.  The  current  flow  and  stored 
energy  in  each  component  vary  with  different  applied  frequencies.  Likewise,  a 
conjunctive-use  system  consists  of  storage  "components"  whose  optimal  control 
decisions  depend  on  the  varying  stream-flow  "signal".  The  optimal  release  and 
water  level  of  each  storage  mechanism  vary  with  different  stream-flow 
autocorrelation. 

Changes  in  the  stream-flow  signal  are  most  significant  when  they  affect 
optimal  decisions  that  ration  and  allocate  water  supplies.  By  incorporating 
autocorrelation  in  the  model  of  a  conjunctive-use  system,  we  can  determine 
efficient  hedging  policies  that  are  more  efficient  and  useful  for  practical 
problems.  This  is  an  important  development  of  this  work,  made  possible  by  my 
concurrent  development  of  the  optimization  tool  presented  in  Chapter  IV. 

3.  Literature  Review  of  Conjunctive-Use  Control  Problems 

Nearly  all  early  conjunctive-use  literature  reviewed  by  Burges  and  Maknoon 
[1975]  assumes  that  one  or  several  parameters  dominate  the  problem  at  hand. 
This  assumption  was  necessary  because  of  optimization  methods  available  at  the 
time  were  limited.  Though  methods  and  computational  abilities  have  advanced 
since,  more  recent  problems  remain  excessively  simplified. 

Among  the  earliest  optimal  control  studies  of  conjunctive  use,  Buras  [1963] 
models  a  steady-state  two-component  surface/groundwater  system  using  a 
coarse  discrete  dynamic  programming  (DDP)  approach.  He  models  the  aquifer 
storage  mechanism  essentially  the  same  as  an  equivalent  reservoir.  Pumping 
costs — though  included — are  linear  with  pumping  rate.  The  problem  does  not 
consider  the  quadratic  effect  of  draw  down  on  cost  or  rate  limits  on  pumping  or 
recharge. 

Buras  [1972]  continues  this  early  work,  proposing  but  not  solving  a 
stochastic  four-dimensional  sequential  decision  problem  that  includes  a 
reservoir,  aquifer  storage,  salinity,  and  prior  month's  stream  flow.  To  reduce 
memory  requirements,  he  recommends  a  polynomial  return  function  to 
approximate  the  minimum  expected  cost  of  future  system  operations.  A 
heuristic  rule  specifies  control  of  groundwater  and  restricts  it  to  a  backup  role 
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used  only  after  the  reservoir  is  empty  or  full.  A  constant  polynomial  function  of 
pumping  rate  approximates  the  changing  pumping  costs  resulting  from 
groundwater  draw  down. 

Gal  [1979]  focuses  on  the  autocorrelation  of  inputs  by  proposing  a 
conjunctive-use  model  for  Lake  Kinneret  and  two  aquifers  that  make  up  a  major 
part  of  Israel's  system.  Gal  incorporates  two  prior  stream-flow  measurements 
using  a  second-order  auto-regressive  model.  Because  this  model  poses  a  difficult 
optimal  control  problem  with  five  dimensions,  he  simplifies  the  model  and 
solves  an  optimal  control  problem  with  three-dimensions  that  does  not  include 
the  aquifers.  Also,  Gal's  proposed  problem  does  not  include  constraints  on 
pumping  rate  and  assumes  linear  rationing  and  pumping  costs. 

Lettenmaier  and  Burges  [1979]  investigate  the  ranges  of  storage  capacity, 
groundwater  pumping  and  recharge  parameters  for  which  conjunctive-use 
appears  feasible.  Stochastic  stream  flows  are  incorporated  by  a  Monte-Carlo 
analysis  of  synthetic  stream-flow  series.  A  heuristic  rule  is  used  to  specify 
control  of  groundwater,  restricting  it  to  a  backup  role.  Lettenmaier  and  Burges 
convert  the  effective  storage  benefit  obtained  by  aquifer  storage  to  an  equivalent 
surface  reservoir,  and  observe  that  effective  storage  decreases  as  pumping  and 
recharge  capacities  are  reduced.  However,  they  note  that  the  effect  is  not 
substantial  until  these  capacities  are  reduced  below  about  one-half  the  annual 
system  demand. 

Bogle  and  O'Sullivan  [1979]  study  the  effect  of  hedging  in  a  system 
consisting  of  a  reservoir,  uncertain  stream  flow,  and  an  alternative  source.  The 
alternative  source  provides  water  at  a  fixed  price  per  unit  volume  which  could 
approximate  a  groundwater  source  under  certain  conditions.  Though  the 
alternative  source  can  provide  water  with  certainty,  it  alone  cannot  meet  the 
whole  demand,  so  shortfall  is  possible.  The  authors  specify  the  supply  rule  for 
the  alternate  source  explicitly  as  the  difference  between  demand  and  release,  up 
to  the  alternative's  maximum  capacity.  Costs  arise  from  use  of  the  alternate 
source  and  penalty  cost  for  low  storage  and  shortfall.  Hedging  occurs  in 
decisions  that  attempt  to  maintain  storage  levels,  though  the  approach  is 
heuristic. 

Existing  conjunctive-use  solutions  that  consider  stochastic  inputs 
generally  rely  on  heuristic  control  policies.  Heuristic  control  policies  used  by 
Buras  [1972],  Lettenmaier  and  Burges  [1979],  and  Bogle  and  O'Sullivan  [1979]  fail  to 
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consider  the  possible  benefits  of  more  flexible  operating  guidelines  that  might 
anticipate  and  mitigate  future  rationing  costs. 

Willis  and  Yeh  [1987]  presented  two  complex  conjunctive-use  models  that 
incorporate  many  realistic  problem  characteristics.  Their  first  model  is  a  two- 
component  surface/ groundwater  system  that  uses  a  response  matrix  approach 
to  incorporate  groundwater  heads,  allowing  the  use  of  quadratic  pumping  costs 
and  bounds  on  heads.  Their  second  model  is  a  three-period  linear  problem 
containing  49  constraints,  24  decision  variables,  and  30  state  variables  using  the 
simplex  method  and  based  on  the  Arkansas  River  Valley  in  Colorado.  Both 
problems,  however,  avoid  stochastic  inputs. 

The  conjunctive-use  problem  encompasses  two  areas  of  research  in  the 
modeling  and  management  of  hydrological  systems:  surface  reservoir  control 
and  aquifer  management.  Uncertainty  from  stream-flow  inputs  presents  the 
largest  challenge  in  determining  control  decisions  for  reservoir  control  problems. 
In  contrast,  parameter  uncertainty  is  the  principle  focus  of  many  aquifer 
management  problems.  Though  input  uncertainty  and  parameter  uncertainty 
both  affect  optimal  control  decisions,  they  are  very  different.  Input  uncertainty  is 
dynamic  and  makes  it  impossible  to  identify  a  single  optimal  control  schedule  in 
advance.  In  contrast,  parameter  uncertainty  is  not  dynamic:  once  quantified, 
this  uncertainty  does  not  evolve  with  the  system.  For  a  problem  that  considers 
only  parameter  uncertainty  and  not  input  uncertainty,  we  can  still  identify  a 
single  optimal  control  schedule. 

Researchers  in  reservoir  control  and  aquifer  management  have  used  two 
different  approaches  in  solving  problems.  On  one  hand,  reservoir  control 
problems  reduce  all  systems  to  simple  lumped  parameter  models  that  do  not 
consider  smaller  scale  variability  in  the  phenomena  described.  On  the  other 
hand,  aquifer  management  problems  generally  ignore  input  uncertainty  and  the 
effect  this  uncertainty  has  on  the  control  of  dynamic  systems. 

These  two  approaches  have  developed  because  of  limitations  on  available 
optimization  techniques.  By  ignoring  input  uncertainty,  we  can  solve  problems 
by  a  variety  of  optimization  techniques  and  can  model  complex  systems  using 
many  variables  [Willis  and  Yeh,  1987].  By  including  uncertainty,  we  generally  are 
restricted  to  models  with  few  variables  and  to  the  optimization  technique  of 
discrete  dynamic  programming  (DDP:  not  the  same  as  differential  dynamic 
programming)  that  quickly  becomes  computationally  infeasible  with  additional 
variables  [Yakowitz,  1982;  Yeh,  1985].  In  general,  aquifer  management  problems 
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ignore  input  uncertainty  [Gorelick,  1983;  Casola  et  al,  1986;  Gorelick,  1987;  Jones  et 
al,  1987;  Wagner  and  Gorelick ,  1987, 1989;  Willis  and  Yeh,  1987;  Yazicigil  and 
Rasheeduddin,  1987;  Georgakakos  and  Vlatsa,  1991;  Wagner  et  al,  1992;  Tiedeman  and 
Gorelick,  1993]  and  most  reservoir  control  problems  rely  on  simplified  models 
[Yakowitz,  1982;  Yeh,  1985;  Saad  and  Turgeon,  1988;  Kelman  et  al, 1990;  Johnson  et  al, 
1991;  Karamouz  and  Vasiliadis,  1992;  Saad  et  al,  1992].  Maddock  [1974]  developed 
operating  rules  for  a  stream-aquifer  system  that  incorporates  uncertain  demand, 
though  the  problem  solves  only  a  single  optimal  control  schedule  based  on  an 
equivalent  deterministic  formulation. 

In  this  work,  I  use  an  optimization  method  that  is  more  flexible  than 
methods  available  for  past  work.  This  allows  me  to  solve  a  conjunctive-use 
problem  characterized  by  stochastic  inputs  with  greater  accuracy  and  less 
simplification  than  required  by  past  studies.  Chapter  IV  develops  and  discusses 
the  optimization  method  of  Gradient  Dynamic  Programming. 

Though  a  basic  conjunctive-use  system  is  relatively  simple,  prior  efforts 
have  not  addressed  two  important  characteristics  of  these  systems.  The  first  is 
the  effect  that  autocorrelation  can  have  on  optimal  control.  This  research  is  the 
first  to  consider  this  effect  without  heuristic  constraints  on  system  control.  This 
allows  us  to  observe  control  decisions  that  hedge  and  take  advantage  of  the 
distinct  characteristics  of  reservoir  and  aquifer  storage.  The  second  is  the  effect 
that  distributed  characteristics  of  groundwater  systems  can  have  on  pumping 
costs.  Though  not  addressed  in  this  thesis,  I  hope  to  address  this  second 
characteristic  in  later  work  by  linking  distributed  parameter  models  of  aquifer 
management  with  the  discrete  parameter  models  of  reservoir  control,  perhaps  in 
an  approach  similar  to  Willis  and  Yeh,  1987. 
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B.  OPTIMAL  CONTROL 


1.  Concept 

There  are  a  variety  of  criteria  that  can  be  used  to  evaluate  the  performance 
of  control  decisions  applied  to  the  management  of  a  conjunctive-use  system. 
These  criteria  include  maximum  revenue,  minimum  risk,  or  minimum  deviation 
from  a  standard.  In  this  problem,  performance  is  measured  by  the  cost  of  control 
decisions;  and  optimal  control  decisions  are  those  that  minimize  expected  cost. 
Cost  is  a  convenient  and  common  basis  for  valuing  an  objective  and  is  especially 
convenient  when  combining  different  objectives  in  a  multi-objective  problem. 

For  simplicity,  I  assume  that  all  costs  are  known  and  quantified.  Chapter  III 
discusses  the  formulation  of  an  objective  function  when  costs  are  ambiguous  or 
difficult  to  quantify. 

The  manager  of  a  water  system  must  periodically  make  two  types  of 
decisions  that  determine  how  water  is  allocated.  The  first  is  how  to  allocate 
water  in  "space" — i.e.,  between  the  storage  mechanisms  and  uses — while  obeying 
system  constraints  and  anticipating  the  future  stresses  on  the  system.  The  second 
is  how  to  allocate  water  in  time;  when  and  how  much  to  allocate  to  demands  to 
best  achieve  benefits  or  avoid  costs. 

Decisions  that  allocate  water  in  space  and  in  time  are  both  affected  by 
uncertainty.  A  manager’s  need  to  anticipate  the  system's  response  without 
certain  knowledge  of  future  hydrological  conditions  leads  to  the  decision  to 
hedge.  For  a  water  supply  system,  a  manager  might  hedge  by  rationing  supplies 
to  reduce  the  chance  of  more  severe  shortages  in  the  future.  For  a  flood  control 
problem,  a  manager  might  hedge  by  allowing  limited  short  term  flood  damage  to 
reduce  the  chance  of  more  severe  flooding  later.  These  are  examples  of  how  a 
manager  may  hedge  allocation  decisions  in  time.  Though  incurring  extra  costs 
for  recharge  and  pumping,  a  manager  might  hedge  by  storing  water  in  the 
ground  to  allow  greater  opportunity  to  retain  flood-water  in  a  reservoir, 
simultaneously  maximizing  opportunity  to  use  flood  water  and  minimize  flood 
risks.  This  is  an  example  of  how  a  manager  might  hedge  allocation  decisions  in 
space. 

A  manager  makes  periodic  "real-time"  control  decisions  based  on 
currently  available  information  about  the  system.  System  operators  apply  these 
decisions  for  a  discrete  period,  or  a  "stage",  of  operations  until  a  manager  verifies 
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or  updates  control  decisions  based  on  changes  in  system  information.  These 
decisions  are  based  on  various  unchanging  characteristics  of  the  system  that 
constrain  control  decisions  (such  as  reservoir  sizes,  aqueduct  capacities,  statutory 
requirements,  contracts,  etc.)  and  on  various  changing  characteristics  that  affect 
how  one  operates  within  the  constraints  (such  as  current  reservoir  levels, 
demands,  stream  flows,  snow-pack  accumulation,  weather  patterns,  time  of  the 
year,  etc.).  Though  both  the  changing  and  unchanging  characteristics  of  the 
system  affect  the  optimum  control  decisions,  only  those  characteristics  that 
change  over  time  produce  changes  in  the  control  decisions.  These  changing 
characteristics  of  the  system  describe  the  state  of  the  system — given  the  state  of 
the  system,  one  knows  everything  that  is  available  to  help  identify  optimal 
control  decisions.  Both  the  state  of  the  system  and  the  control  decisions  can  be 
summarized  by  variables  that,  respectively,  are  called  state  variables  and 
decision  variables.  The  goal  of  the  optimal  control  problem  is  to  determine  the 
optimal  value  of  decision  variables  as  a  function  of  state  variables. 

2.  General  Model 

To  determine  optimal  controls  for  a  real  system,  we  may  describe  the 
system  by  an  abstract  mathematical  model  that  describes  the  characteristics, 
constraints,  and  objectives  of  system  operation.  In  the  work  presented  here  I 
assume  that  a  system  can  be  represented  by  a  mathematical  model  that  lumps 
system  characteristics  into  a  limited  number  of  variables. 

The  following  generic  model  uses  standard  "state-space"  notation 
frequently  used  in  reservoir  management  modeling. 

a.  Variables 

Control  decisions  are  described  by  a  vector  of  m  "decision  variables"  u. 
These  decisions  may  include  how  much  water  to  release  to  meet  needs,  how  to 
allocate  available  supplies  to  available  storage,  and  so  forth.  A  manager  makes 
these  decisions  based  on  current  values  of  changing  characteristics  of  the  system 
defining  the  initial  "state"  of  the  system  and  described  by  a  vector  or  n  "state 
variables"  x.  The  system  evolves  to  state  y  at  the  end  of  a  stage  with  the 
application  of  control  decisions  and  realizations  of  stochastic  inputs.  Stochastic 
inputs  are  represented  by  a  vector  of  r  "stochastic  variables"  w.  They  can 
include  uncertain  stream  flows,  measurement  errors,  or  other  system 
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characteristics  that  prevent  exact  prediction  of  y.  This  evolution  is  modeled  by  a 
transition  function 

y  =  T(x,u,w)  . 

For  simplicity,  parameters  x,  u,  w,  and  y  are  with  reference  to  the  current 
stage  and  thus  are  not  indexed.  A  state  y  at  the  end  of  a  stage  t  is  the  state  x 
at  the  beginning  of  the  next  stage  f+1;  i.e.,  yt  =  xt+i . 

Though  the  actual  values  of  w  are  unknown  until  the  end  of  the  stage,  I 
assume  that  the  multivariate  probability  distribution  of  w  is  known  and 
composed  of  independent  stochastic  variables.  When  w  contains  variables  that 
are  autocorrelated,  such  as  time  series  of  stream  flows,  realized  values  of  w 
from  prior  stages  condition  the  probability  distribution.  These  prior  values  also 
affect  the  expected  cost  of  control  decisions  and  the  optimal  solution  u*(x);  and 
as  important  information  on  the  state  of  the  system,  they  are  included  as  state 
variables.  With  this  extension  of  x,  we  can  specify  the  distribution  of  w  by  an 
unchanging  or  "stationary"  model  w(x)  that  depends  only  on  the  state  vector. 

b.  Constraints 

Without  constraints,  variables  describing  a  system  can  take  on  values  that 
are  not  feasible  for  a  real  system.  Some  constraints  exist  due  to  relationships 
between  variables  such  as  those  defined  above;  some  exist  due  to  physical  or 
practical  limits.  Frequently,  constraints  exist  on  the  potential  states  of  the  system 
and  on  control  decisions.  These  are  bounds  on  state  and  decision  variables: 

Xmin  —  X  ^  Xmax  ,  ymin  —  y  —  ymax  >  2nd  Uniin  —  ^  —  Umax  • 

Feasibility  of  initial  state  x  is  ensured  by  choosing  only  feasible  values;  thus 
constraints  on  x  do  not  enter  into  the  optimization  method  discussed  below. 
Other  constraints  may  exist  to  define  variables  needed  to  model  a  system  but  are 
represented  by  equality  constraints  in  the  transition  equation  or  directly  in  the 
objective  function. 

Constraints  must  be  respected  by  control  decisions.  The  standard  form  of 
inequality  constraints  is 


g(x,u)  <  d(x) 
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where  g  and  d  are  scalar  functions.  Using  the  transition  equation 
y  =  T(x,u,w)  ,  we  can  formulate  constraints  that  do  not  depend  on  y .  I  assume 
drat  decisions  must  be  feasible  regardless  of  any  realized  value  of  stochastic 
input  w;  thus,  with  only  bounds  as  constraints  on  variables,  we  can  more  simply 
represent  them  by 

gTu  <  d(\) 

where  g  is  a  vector  of  scalar  coefficients.  Appendix  A  provides  further 
discussion  of  the  development  of  constraints  for  reservoir  control  problems. 

c.  Objective  Function 

The  measure  of  performance  of  the  decisions  is  an  objective  function  f 
that  measures  expected  cost  of  control  decisions.  Costs  might  include  losses 
from  hedging  (e.g.,  rationing  and  minor  short-term  flood  or  environmental 
damage),  operational  costs /benefits  (e.g.,  allocation  expenses  such  as  from 
pumping  and  hydropower  generation)  or  a  variety  of  qualities  that  produce 
extra-market  costs  or  benefits  (e.g.,  water  quality,  recreation,  and  wildlife 
habitat).  The  objective  function  can  be  dependent  on  control  decisions,  state  of 
the  system,  and  uncertain  inputs.  Given  the  "optimal  decisions"  u*  =  u*(x)  that 
minimize  costs  given  a  state  x,  the  optimum  value  of  the  objective  function  is 

/*(x)  =  min„{/(x,u,w,y)}  =/(x,u*,w,y) 

For  optimal  control  of  reservoir  systems,  Foufoula-Georgiou  and  Kitanidis 
[1988]  propose  a  general  objective  function  composed  of  two  costs:  those 
incurred  in  the  current  stage  and  those  expected  in  following  stages.  Current 
cost  results  directly  from  current  control  decisions  and  may  also  depend  on  the 
state  of  the  system.  Expected  cost  also  results  indirectly  from  current  control 
decisions  because  of  the  effect  that  current  decisions  have  on  limiting  or 
changing  future  decisions  and  states  of  the  system.  The  objective  function  is 
represented  by 


/=  Ew{C(x,u,w)  +  F(y)} 

where  C  is  the  current  cost  function  and  F  is  the  expected  cost  function  (also 
known  as  the  "cost-to-go"  from  state  y).  Because  w  is  uncertain,  y,  F,  and  C 
are  also  uncertain.  The  objective  function  is  the  expected  cost  calculated  over  the 
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distribution  of  w  and  represented  by  the  operator  Ew{}.  This  is  estimated  by 
the  probability  weighted  sum  of  costs  calculated  from  discrete  realizations  w(*) , 
thus 


/  =  X  {p*[C(x,u,w(*))  +  F(y(jt))]} 
k 

for  weights  pk ,  where  y(*>  =  T(x,u,ww)  . 

This  function  allows  inclusion  of  uncertain  costs  in  the  current  stage; 
however,  this  can  be  simplified  to 

/  =  C(x,u)  +  2  {p*F(y(*))) 

k 

to  consider  only  deterministic  current  costs  by  choosing  stages  that  are 
sufficiently  short.  Foufoula-Georgiou  and  Kitanidis  [1988]  also  use  this  practical 
simplification. 
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C.  PROBLEM  DESCRIPTION 


The  sample  problem  models  and  optimizes  control  decisions  for  a 
conjunctive-use  system.  The  problem  contains  no  heuristic  control  rules  that 
restrict  the  management  of  groundwater,  and  thus  the  optimal  control  solution 
has  a  greater  set  of  feasible  decisions  from  which  to  find  the  optimal.  This  allows 
us  to  observe  the  effect  that  the  different  capabilities  and  constraints  of  reservoir 
and  aquifer  storage  have  on  control  decisions,  especially  with  respect  to  hedging 
of  supplied  water  and  allocation  of  stored  water.  In  addition  we  can  observe 
how  the  persistence  of  stream  flows  affect  hedging  and  allocation. 

The  sample  problem  uses  a  generic  conjunctive  use  model  similar  to  that 
discussed  by  Buras  [1963, 1972]  and  summarized  by  Yakowitz  [1982].  It  is  based 
on  the  EBMUD  system  that  supplies  water  to  communities  along  the  eastern 
edge  of  the  central  and  southern  San  Francisco  Bay.  Figure  II-l  illustrates  the 
EBMUD  system  and  conjunctive  use  model. 

1.  System  Description 

The  EBMUD  system  serves  the  residential  needs  of  approximately  1.2 
million  people,  as  well  as  the  industrial,  commercial,  and  institutional  needs  in 
the  East  Bay  region  of  the  San  Francisco  Bay  area.  Notable  cities  in  the  area 
include  Oakland  and  Berkeley.  About  95  percent  of  its  water  supply  is  from  the 
Mokelumne  River’s  575-square  mile  watershed  on  the  western  slope  of  the  Sierra 
Nevada  Mountains. 

Source 

Supplies  from  stream  flow  are  seasonal  and  uncertain,  with  flow  into 
EBMUD’s  reservoir  system  at  about  720  thousand  acre-feet  (TAF)  annually 
(based  on  United  States  Geological  Survey  records  since  1928  for  the  Mokelumne 
Hill  Gage).  Upstream  of  EBMUD's  system  there  are  small  diversions  of  about  10 
TAF  annually  by  rural  Amador  and  Calaveras  Counties.  Annual  flow  has  varied 
from  a  high  of  1,595  TAF  in  1983  to  a  low  of  130  TAF  in  1977.  Average  monthly 
flows  vary  from  a  high  of  over  100  TAF  in  May  to  a  low  of  about  30  TAF  through 
the  fall  months.  However,  a  typical  year  will  generally  have  a  month  of  high 
flow  approaching  200  TAF  and  a  month  of  low  flow  approaching  10  TAF;  and,  at 
times,  the  natural  flow  of  the  river  has  ceased. 
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Knowledge  of  monthly  flow  distributions,  autocorrelation,  and  snow 
accumulation  provides  valuable  information  for  improving  predictions  of  flow. 
Intra-annual  stream  flows  are  highly  autocorrelated.  Also,  much  of  the  late 
spring  and  summer  runoff  is  a  result  of  melting  snowpack,  and  measurements  of 
snowpack  are  available  throughout  the  winter.  Therefore,  consideration  of  the 
season  and  measurements  of  snowpack  and  prior  stream  flows  all  contribute  to 
improvements  in  stream-flow  prediction. 

Existing  System 

EBMUD  manages  two  reservoirs  having  a  combined  capacity  of  641 
thousand  acre-feet  (TAF)  on  the  Mokelumne  River.  Up  to  200  TAF  of  this  storage 
is  reserved  for  flood  control  under  agreements  with  the  U.  S.  Army  Corps  of 
Engineers,  and  an  additional  21  TAF  is  dead  storage.  Besides  water  supply  and 
flood  control  functions,  the  reservoirs  also  have  a  combined  hydropower 
generating  capacity  of  39  megawatts. 

An  82  mile  aqueduct  transports  water  to  the  service  area  by  for  use  or 
storage  in  five  terminal  reservoirs.  The  terminal  reservoirs  provide  an  additional 
155  TAF  of  storage  capacity  (with  17  TAF  of  dead  storage).  The  aqueduct  has  a 
delivery  capacity  of  200  million  gallons  per  day  (MGD)  by  gravity  flow  or  325 
MGD  by  pumping.  This  maximum  delivery  capacity  also  coincides  with 
EBMUD’s  permits  to  divert  up  to  325  MGD  (364  TAF  per  year)  for  use  in  its 
service  district. 

Proposed  Aquifer  Storage 

Both  the  River  and  Aqueduct  flow  west  from  the  Sierras  and  across  the 
Central  Valley  of  California  before  meeting  the  San  Francisco  Bay  Delta.  This 
area  contains  large,  extensive  aquifers. 

Underlying  this  area,  the  fresh-water  bearing  Victor,  Laguna,  and  Mehrten 
formations  generally  consist  of  thick  Quaternary  and  late  Tertiary  sand  and 
gravel  deposits.  Under  these  are  brackish  and  saline  marine  sediments.  Though 
locally  confined  conditions  exist,  the  upper  formations  generally  form  an 
unconfined  aquifer.  Thicknesses  of  these  formations  total  a  few  hundred  feet  in 
the  east  increasing  to  almost  two  thousand  feet  near  the  Delta.  Well  capacities 
have  been  reported  in  the  range  of  500  to  1500  gallons  per  minute  (gpm),  and 
studies  indicate  specific  capacities  in  the  range  of  35  to  60  gpm/ft  and 
transmissivities  of  60,000  to  80,000  gal / day/ ft  [D WR,  1967]. 
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Due  to  pumping  for  local  agricultural  and  municipal  needs,  the  water 
table  is  depressed  by  an  average  50  to  100  feet  below  virgin  levels.  Currently  the 
water  table  averages  about  100  feet  below  land  surface  and  is  below  sea  level  in 
many  locations.  This  has  resulted  in  salt-water  intrusion  near  the  Delta  and  has 
prompted  local  agencies  to  seek  remedies  that  reduce  over-pumping  in  the  area. 
However,  as  a  result  of  this  over  pumping,  there  is  a  great  amount  of  available 
space  for  aquifer  storage. 

Additional  aquifer  storage  opportunities  exist  in  the  EBMUD  service 
district  and  other  areas  adjoining  the  District.  Indeed,  neighboring  Alameda 
County  Water  District  already  has  a  local  program  for  aquifer  replenishment, 
though  its  efforts  are  aimed  at  restoring  groundwater  levels  and  preventing  salt 
water  intrusion  rather  than  at  managing  the  aquifer  as  a  storage  mechanism  [Sax 
et  al,  1991].  Nevertheless,  EBMUD  has  identified  the  lower  Mokelumne  region 
in  the  Central  Valley  as  the  best  site  based  on  storage  capacity,  well  capacities, 
water  quality,  institutional  constraints,  and  anticipated  costs  [EBMUD,  1992, 
Appendix  D3]. 

2.  Optimal  Control  Problem 

The  conjunctive-use  model  is  a  simplified  representation  of  the  EBMUD 
system  with  an  added  aquifer-storage  component.  The  model  does  not  solve 
specific  management  or  planning  problems  for  EBMUD,  but  it  does  illustrate 
some  of  the  general  characteristics  of  optimal  control  of  a  conjunctive-use  system 
and  points  to  differences  between  reservoir  and  aquifer  storage. 

System  Model 

I  model  the  system  with  three  state  variables.  These  represent  reservoir 
storage  level  (xi),  aquifer  storage  level  fo),  and  prior  month's  stream  flow  (X3). 
Control  decisions  specify  the  amount  of  water  supplied  («i)  to  meet  demands 
and  the  change  in  groundwater  storage  (M2)  resulting  from  recharging  and 
pumping  the  aquifer.  These  decisions  are  applied  over  stages  that  are  one  month 
long.  During  each  stage,  the  state  evolves  from  beginning  values,  x,  to  the 
ending  values,  y,  based  on  control  decisions  and  realized  value  of  uncertain 
stream  flow  (wi ). 

This  effort  is  the  first  to  illustrate  the  effect  that  the  different  capabilities 
and  constraints  of  reservoir  and  aquifer  storage  have  on  hedging  and  allocation, 
without  the  use  of  heuristic  rules  controlling  the  use  of  groundwater.  In 
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addition,  this  effort  includes  the  effect  of  autocorrelation  and  illustrates  the 
influence  that  autocorrelation  can  have  on  hedging  and  allocation. 

Nevertheless,  a  practical  and  complete  model  of  EBMUD's  potential 
conjunctive-use  system  is  more  complex,  with  more  state  and  decision  variables, 
than  presented  here.  A  more  complete  model  could  require  six  or  more  state 
variables  (Table  1).  Here  I  use  only  three;  I  limit  my  current  effort  primarily  to 
illustrate  more  easily  a  general  conjunctive-use  modeling  and  solution  approach 
and  to  determine  some  general  properties  of  conjunctive-use  solutions. 

Transfer  Equation 

Reservoir  levels  change  linearly  with  stream  flow  and  with  control 
decisions  that  supply  and  allocate  water.  Groundwater  level  changes  with  the 
allocation  decision  that  recharges  or  pumps.  The  transfer  equation  incorporates 
these  effects  as  follows 
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Constraints 

I  incorporate  fundamental  differences  between  reservoir  and  aquifer 
storage  by  adding  constraints  on  control  decisions.  As  discussed  above  aquifers 
generally  have  much  greater  storage  capacities  than  reservoirs  but  more  limited 
response  rates  from  filling  and  extraction. 

We  may  be  able  to  design  reservoirs  with  storage  capacities  approaching 
that  of  many  available  aquifers;  we  may  also  be  able  to  design  pumping  and 
recharge  facilities  that  allow  rapid  changes  in  groundwater  levels.  However,  a 
system  with  reservoir  and  aquifer  storage  can  make  effective  use  of  both  storage 
mechanisms  without  building  facilities  needed  to  remove  these  constraints.  By 
efficiently  combining  the  two  storage  mechanisms,  we  can  achieve  a  high  level  of 
water  supply  reliability  with  much  less  expense. 

Table  2  summarizes  the  constraints  that  I  apply  to  the  conjunctive-use 
problem.  These  constraints  account  for  differences  between  reservoir  and  aquifer 
storage  by  specifying  maximum  reservoir  level  and  maximum  change  in  aquifer 
level  during  a  one-month-long  stage.  In  standard  form,  these  constraints  are: 

(a)  Upper  bound  of  reservoir  capacity:  [-l,-l]u  £  400  -  x\  -  w\ 
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(b)  Lower  bound  of  reservoir  capacity: 

(c)  Lower  bound  of  aquifer  storage: 

(d)  Upper  bound  on  recharge: 

(e)  Upper  bound  on  pumping: 

(f)  Upper  bound  on  demand: 

(g)  Lower  bound  on  demand: 


[  1,  l]u  <  Xi  +Wi 
[  0,-1  ]u  <  x2 
[0,  l]u  <  10 
[  0,-1  ]u  <  20 
[  1, 0]u  <  58.3 
[-l,0]u  <  0 


where  numbers  represent  thousands  of  acre-feet  (TAF). 

Demand  is  assumed  to  be  700  TAF  per  year.  This  is  close  to  the  mean 
annual  stream  flow.  Thus,  demand  in  each  stage  is  58.3  TAF  since  control 
decisions  are  updated  monthly. 

The  stochastic  variable  in  these  constraints  presents  a  difficulty.  An 
inequality  constraint  becomes  an  equality  constraint  when  binding,  and  the 
presence  of  a  stochastic  variable  in  a  constraint  means  that  there  is  not  a  unique 
control  solution  that  will  satisfy  the  constraint  for  each  possible  realization  of  the 
stochastic  variable. 

I  use  two  assumptions  to  address  the  difficulty  of  stochastic  variables  in 
the  constraints.  This  first  assumption  is  that  demands  can  be  met  only  by  stored 
water  and  not  by  current  stream  flow.  As  a  result,  the  supply  and  allocation 
decisions  are  more  conservative  than  necessary,  but  this  avoids  the  difficulty  of 
specifying  decisions  that  change  throughout  the  month  to  adjust  to  changing 
stream  flow.  The  second  assumption  is  that  the  discrete  realization  of  stream 
flow  that  first  produces  a  binding  constraint  is  the  only  value  needed  to  specify  a 
stochastic  constraint.  For  example,  stochastic  constraint  (b)  is  represented  by  a 
set  of  constraints  resulting  from  each  discrete  realization  of  W] : 


gTu  <  Xi+  Wli(i)  ,  gTu  <  X\  +  Wi>(2)  /  gTu  ^  *i  +  wi.(3)  /•••■• 

where  g  =  [1, 1]T.  If  wi,(i)  is  the  smallest  discrete  realization,  only  the  first 
constraint  can  be  become  an  equality  constraint;  converting  any  other  constraint 
in  the  set  to  an  equality  constraint  is  infeasible  since  this  would  violate  the  first 
constraint.  Thus,  w1>(i)  is  the  only  value  needed  to  specify  the  stochastic 
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constraint  (b).  If  stream  flow,  w\ ,  can  attain  a  minimum  value  of  zero  with  non¬ 
trivial  probability,  it  is  reasonable  to  replace  stochastic  constraint  (b)  with 


(h)  [1,  l]u  ^  x\  . 

Control  decisions  will  then  be  feasible  for  any  realized  value  of  stream  flow. 

A  problem  with  these  two  assumptions  is  that  control  decisions  can  be 
excessively  constrained.  It  is  even  possible  that  the  constraints  result  in  no 
feasible  decisions.  To  prevent  this,  we  can  choose  a  decision  period  that  is 
sufficiently  short,  and  we  can  limit  selection  of  discrete  values  W(*)  to  reasonable 
values  that  are  contained  in  a  bounded  set.  In  the  conjunctive-use  problem 
presented  here,  any  non-negative  value  of  W(*)  is  possible,  though  the 
probability  of  extremely  large  values  becomes  insignificant  and  does  not  affect 
the  optimal  solution.  Because  the  objective  function  includes  only  demand  for 
water  (and  not,  for  example,  flood  control),  releases  in  excess  of  demand  required 
to  meet  the  reservoir  capacity  constraint  (a)  have  no  effect  on  the  objective  and 
optimal  controls.  If  a  release  is  required,  the  reservoir  is  full  and  thus  the  amount 
released  does  not  affect  future  decisions  or  costs.  Thus,  for  this  problem,  I  drop 
constraint  (a)  and  set  yi  —  400  whenever  the  transition  equation  results  in 
yi  >  400. 

An  alternate  approach  to  handling  stochastic  constraints  that  I  do  not  use 
here  is  to  add  stochastic  slack  variables  s  to  the  transition  equation.  This  results 
in  the  transition  equation 

'looi  r-i  -ii  1 1  r-r 

y  =  T(x,u,w  )=  010  x+  01  u+  0  w  +  0  s  • 
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The  values  of  slack  variables  depend  on  control  decisions  and  realized  stream 
flow,  thus  s(u,w)  .  Though  more  complicated,  the  solution  algorithm  presented 
in  Chapter  IV  could  be  expanded  to  include  these  variables. 

Objective  Function 

The  objective  function  is  the  recursively  derived  function 

/  =  C(u)  +  X  (P*F(y(*))} 
k 

where  C(u)  =  (5  -  4wi/58.3)(58.3  -  «i)  +  0.0001  u22 
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and 


F(M)(y)  =  °  • 


58.3  is  the  monthly  water  demand.  This  objective  function  only  incorporates  the 
cost  of  water  rationing.  It  does  not  include  other  costs  of  flooding  or 
groundwater  pumping,  or  benefits  of  hydropower  generation.  It  does,  however, 
include  a  small  penalty  for  groundwater  pumping  and  recharge  to  make  the 
objective  strictly  convex  with  respect  to  these  decisions.  This  prevents 
unnecessary  pumping  and  recharge. 

Though  not  based  on  a  practical  assessment  of  rationing  costs,  the 
objective  function  contains  many  characteristics  of  a  practical  function. 

Significant  characteristics  of  any  practical  function  are  increasing  total  and 
marginal  costs  as  rationing  increases  (supply  decreases).  This  is  included  in  the 
objective  function. 

This  function  assumes  that  an  unlimited,  fixed-price  water  supply 
alternative  does  not  exist.  As  a  result,  the  function  is  strictly  convex  with  respect 
to  rationing  decisions,  and  optimal  decisions  prefer  long  periods  of  moderate 
rationing  to  short  periods  of  severe  rationing.  Was  I  to  use  a  linear  function 
representing  the  presence  of  a  fixed-price  water  supply  alternative,  rationing 
would  not  occur  beyond  a  level  that  produces  marginal  costs  equal  to  the  fixed 
price.  In  practical  problems,  this  condition  rarely  exists.  However,  even  if  I  was 
to  consider  this  possibility,  the  objective  function  would  still  be  convex,  but  not 
strictly  convex.  Concave  functions  are  rare;  optimal  control  with  a  concave 
objective  function  can  produce  an  unlikely  all-or-nothing  approach  to  allocating 
supplies. 

Use  of  the  above  objective  function  is  also  justified  since  it  incorporates 
rationing  costs,  the  most  significant  component  of  operating  cost  for  many 
systems.  Most  other  costs  are  less  significant  when  compared  to  rationing  costs. 
For  a  conjunctive-use  system  with  wells  pumping  at  1000  gpm  and  Central 
Valley  aquifer  conditions  discussed  above  in  Section  C.l,  pumping  costs  are 
about  $10  to  $20  per  acre-foot  based  on  just  the  cost  of  electricity  [Georgakakos  and 
Vlatsa,  1991].  For  the  EBMUD  system,  the  historical  value  of  hydropower 
generated  by  the  EBMUD  system  is  about  $10  per  acre-foot  [EBMUD,  1992,  Vol. 
1].  In  contrast,  when  water  is  scarce,  the  marginal  value  of  municipal  water 
supplies  can  be  much  greater.  California  sold  water  from  its  1992  Water  Bank  for 
$175  an  acre-foot  [McClurg,  1992b].  The  California  communities  of  Santa 
Barbara,  Goleta  and  Montecito,  and  others  have  developed  desalination  plants  to 
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provide  backup  water  supplies  that  cost  $1500  to  $2000  an  acre-foot.  Even  when 
water  is  not  scarce,  wholesale  costs  south  of  the  Delta  range  from  $44  in  San 
Joaquin  Valley  to  $237  in  southern  California  [McClurg,  1992a]. 

3.  Stochastic  Model  of  Stream-Flow  Time  Series 

Figure  II-2  displays  the  monthly  stream-flow  hydrograph  of  the 
Mokelumne  River.  I  model  the  probability  distributions  of  monthly  stream  flows 
from  this  time  series  by  using  a  reasonable  three-parameter  model.  The  model 
also  incorporates  an  additional  parameter  for  each  month  to  account  for  stream- 
flow  autocorrelation  with  the  prior  month.  For  a  practical  model  of  the  EBMUD 
system,  a  stream-flow  model  could  also  include  measurements  of  snowpack 
accumulation  during  the  winter  and  spring;  the  approach  used  to  incorporate 
this  additional  information  would  be  similar  to  that  used  to  incorporate 
autocorrelation. 

Model  of  Distribution 

The  probability  distribution  for  each  months  stream  flows  are  modeled  by 
a  seasonal  3-parameter,  log-normal  model.  Chapter  III  discusses  how  seasonality 
is  incorporated  in  the  optimal  control  problem.  Parameter  Xx  shifts  stream-flow 
values  to  fit  more  closely  a  log-normal  distribution  modeled  by  log-mean  fix  and 
log-standard  deviation  cv  thus 

wr  =  x-t  +  exp(£<xT  +  fix) 

for  each  month  x.  Stochastic  variable  £  has  a  standard-normal  distribution. 
These  parameters  are  estimated  from  the  historical  data  as  discussed  below. 
Examples  of  probability  quantiles  for  monthly  flow  are  displayed  in  Figure  II-3. 

Model  of  Autocorrelation 

Autocorrelation  is  first-order  auto-regressive  (Figure  II-4).  Correlation 
coefficients  are  calculated  after  stream-flows  have  been  transformed  because  this 
produces  coefficients  that  are  more  significant  than  those  calculated  prior  to 
transformation.  Thus,  the  dependence  of  wT  on  wT-\  (where  wz.\  =  x\z)  is 
modeled  by  the  relationship 

(Or  =  Pt  +  Pr^-ia t-i  ~  Pr-l)  +  e°V  »  1  "  P*2 

°T-1 
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where  ah  =  ln(wT-^T)  and  cth-i  -  ln(wx.i  -  Xr-i)  •  Parameter  pT  is  the  sample 
correlation  coefficient  for  month  t. 


Parameter  Estimation 

The  correlation  coefficient  is  estimated  from  historical  data  by 

Cov(ctXr,  GVj) 

Pt  =  VVarC^)  Var(GvJ  ' 

Stedinger  [1981]  discusses  a  variety  of  additional  definitions  for  the  correlation 
coefficient  that  also  may  be  appropriate.  Parameter  values  and  <Xj  are  the 
unbiased  sample  estimates  from  the  transformed  historical  data.  The  third  "shift" 
parameter,  Xi>  °f  the  3-parameter  log-normal  distribution  is  estimated  by 

_  wt(q)  wt(  1  -q)  -  Wf(0.5)2 
T  -  vvj(^)  +  W[(\-q)  -  2wt(0.5 ) 

as  recommended  by  Stedinger  [1980],  where  parameter  wt(q)  is  the  sample  value 
for  the  <7*th  quantile.  Stedinger  recommends  choosing  wt(q)  and  wt(\-q)  to  be 
the  largest  and  smallest  values  observed  for  month  t.  Table  3  summarizes  the 
parameters  calculated  for  each  month. 


Considerations 

An  advantage  of  the  three-parameter  log-normal  model  is  that,  for  the 
example  problem,  dw/dx  becomes 


0, 


0, 


Wl-Xr  Or 
x3  '  Xr-\  ®V-1 . 


The  three-parameter  log-normal  model  also  presents  some  possible 
complications.  These  include 

1.  Values  of  a  historical  data  set  that  result  in  wT-Xr  <  8  where  5  is  a  small 
value:  Using  the  recommendation  of  Stedinger  [1980]  ensures  that  the  estimate  of 
Xt  will  be  smaller  than  the  smallest  observed  value  of  wT.  However,  the 
discretized  state  space  may  include  values,  w*iT,  that  could  violate  this  condition. 
Thus,  I  do  not  allow  values  of  Xt  greater  than  -1 .  The  Mokelumne  River  flow 
series  generally  is  slightly  less  skewed  than  a  fitted  log-normal  distribution;  as  a 
result,  the  samples  x-t  are  negative  for  all  months  but  July  (Table  3). 
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2.  Historical  data  that  are  not  positively  skewed:  A  three-parameter  log-normal 
model  cannot  produce  a  negative  skew;  however,  without  constraints  the  above 
parameter  estimates  will  try  to  fit  a  negative  skew.  To  prevent  infeasible 
parameters,  a  probability  distribution  model  must  satisfy  the  restriction 

wt(q)  +  w,(l-q)  -  2w,(0.5)  >  0  [Stedinger,  1980].  As  skew  approaches  zero,  %% 
approaches  (-<*>).  For  the  Mokelumne  River  flow  series,  all  monthly 
distributions  have  a  positive  skew  except  for  September  and  October.  Xr  f°r 
these  months  is  set  to  -600  to  more  closely  approximate  a  skew  of  zero  (Table  3). 

3.  Synthetic  stream-flow  values  that  are  negative:  When  this  occurs,  the 
generated  flow  value  and  dwi/dx$  are  set  to  zero.  Since  Xt  -  -1,  generated  flow 
values  may  be  shifted  to  less  than  zero. 

4.  Synthetic  and  real  stream-flow  values  that  are  greater  than  the  largest  value  of 
the  discretized  range  of  the  distribution:  To  discretize  the  state-space  of  the 
system,  a  bounded  range  of  flows  must  be  specified;  however,  log-normal 
distributions  are  unbounded.  Fortunately,  optimal  decisions  applied  to  water 
management  problems  are  not  sensitive  to  sufficiently  large  stream-flows  for 
which  the  probability  approaches  zero.  This  is  especially  true  for  problems  with 
only  a  water  supply  objective — optimal  decisions  do  not  change  when  stream 
flows  are  greater  than  the  sum  of  storage  capacity  and  demand  since  the  system 
cannot  make  use  of  flows  beyond  this  sum.  Thus,  when  generated  flows  greater 
than  this  sum  occur,  I  reduce  their  value  to  this  sum  and  can  set  gradients  of  the 
objective  function  V*/-  to  zero. 
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Table  1.  State  Variables  Needed  by  Models 

Model  of  Real  System 

Current  Model 

Facilities:  Reservoirs  (2-3) 

Reservoir  (1) 

Groundwater  Storage  (min  2) 

Groundwater  Storage  (1) 

Stochastic  Inputs:  Prior  Stream  Flow  (1) 

Snow-Pack  Measurements  (1) 

Prior  Stream  Flow  (1) 

Total  min  6 

3 

Table  2.  Constraints  on  State  and  Decision  Variables 
(1000’s  of  acre-feet) 

Minimum 

Maximum 

Surface  Reservoir 

Storage 

Change  per  month 

0.0 

unlimited 

400.0 

unlimited 

Sub-Surface  Reservoir 
Storage 

Change  per  month 

0.0 

-20.0 

unlimited 

10.0 

Release  for  supply 

Supply  per  month 

0.0 

58.3 
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400 


Historical  Monthly  Stream  Flow 
(1000's  acre-feet) 


Fitted  Monthly  Flow  Quantiles 


Fig.  II-3.  Quantiles  of  the  fitted  3-parameter  log-normal  distribution.  Displayed 
quantiles  are  stream-flow  values  resulting  from  the  back-transformation  of  the 
first  and  second  standard  deviates  of  the  normalized  residuals,  where 

wife)  =  Xr  +  exp(//x+  eaT) 
for  each  month,  x,  and  deviate: 

e=-2,  -1,0,  1,2 

and  resulting  quantiles: 

P{wT <*>,(£)}  =  0.054,  0.242,  0.5,  0.758,  0.946 
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Autocorrelation  and  Partial  Autocorrelation 
of  Monthly  Stream  Flows 


Fig.  II-4.  Lag-one  and  lag-two  autocorrelation  coefficients,  p^l)  and  pt(2),  and 

lag-two  partial  autocorrelation  coefficient,  n^l),  for  each  month  r.  Partial 
autocorrelation  is  the  correlation  that  remains  after  accounting  for  shorter  lags. 
Thus  nfi)  =  p^l),  and  afi)  =  [pfi)  -p^DPf-jflMl  -p^DPf-jtt)]- 
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Table  3.  Stream-flow  Parameters 


Parameters  that  transform  stream  flows,  w,  to  a  standard-normal 
distribution  using  the  3-parameter  lognormal  model: 

e  =  [ln(w-^)  -  n]/o 


shift 

parameter 

normalized 

mean 

std  dev. 

lag-1  a.c. 

Month 

X 

P 

<7 

P 

Oct. 

-600  * ** 

6.448 

.018 

.928 

Nov. 

-6.806 

3.617 

.527 

.812 

Dec. 

-1.106 

3.627 

.736 

.776 

Jan. 

-3.512 

3.824 

.698 

.807 

Feb. 

-3.029 

3.875 

.723 

.763 

Mar. 

-10.101 

4.227 

.573 

.821 

Apr. 

-24.363 

4.569 

.477 

.824 

May 

-90.159 

5.271 

.377 

.896 

Jun. 

-35.207 

4.824 

.537 

.821 

Jul. 

** 

3.573 

.612 

.700 

Aug. 

-89.184 

4.804 

.094 

.917 

Sep. 

-600  * 

6.447 

.016 

.935 

*  distribution  has  negative  skew,  %  value  at  -600 

**  value  adjusted  to  -1 


D.  SOLVING  THE  CONJUNCTIVE-USE  PROBLEM 
1.  Method 

The  optimization  method  used  to  make  iterative  improvements  on  the 
solution  is  a  well-established  optimization  technique  known  as  the  method  of 
successive  approximations  [Esmaeil-Beik and  Yu,  1984].  This  method  solves  the 
recursive  relationship 


F(M)(y(M)>  =  Ao(x«) 

where  =  min„{C(x,u)  +  £  (P*F(o(y)J  ) 

k 

backwards  in  time  from  a  specified  expected  cost  function  F^iy).  This  cost 
function  specifies  expected  costs  beyond  the  end  of  a  planning  horizon,  N  stages 
long.  Note  that  y(_i  =  x, . 

If  we  have  sufficient  experience  with  solving  a  particular  class  of  problems 
and  know  the  solutions  present  some  simple  characteristic  functional  form,  we 
can  describe  F  as  a  function  of  y  without  using  the  piece-wise  interpolation 
method  used  here.  The  problem  solving  would  then  consist  of  a  few  fitted 
parameters.  Often,  however,  we  cannot  specify  a  simple  characteristic  functional 
form  unless  we  are  willing  to  make  some  very  broad  assumptions. 

Instead,  we  can  specify  the  value  of  F  at  a  number  of  discrete  states  y(t) 
sufficient  to  allow  us  to  estimate  the  continuous  function  by  interpolation.  Using 
this  piece-wise  continuous  expected  cost  function,  we  can  solve  the  recursive 
relationship  presented  above  by  iteratively  improving  control  decisions  until  the 
optimum  is  identified  (Figure  II-5). 

2.  Interpretation 

Suppose  that  we  are  extremely  short-sighted  in  our  management  of  a 
system,  and  we  chose  the  goal  of  minimizing  cost  for  just  the  current  stage.  By 
ignoring  future  costs,  we  are— in  effect— saying  that  our  expected  cost  function  is 
zero  (or  constant)  for  all  final  states  of  the  system:  F^( y)  =  0.  Optimal  decisions 
u*  are  those  that  give 

Av)(x)  =  min„{C(x,u))  . 
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If  we  determine  the  optimal  decisions  and  cost  for  discrete  initial  states  x©,  we 
can  specify  a  piece-wise  continuous  optimum  control  function  u*(x)  and 
optimal  cost  function  /*(at)(x)  for  all  states  x. 

Now  suppose  that  we  want  to  be  a  little  less  short  sighted,  and  we  decide 
to  minimize  the  total  costs  incurred  over  both  the  current  stage  and  a  following 
stage.  In  the  first  stage  of  this  two  stage  problem,  optimal  decisions  are  those 
that  give 

J*(N- i)(x)  =  minu{C(x,u)  +  X  iPkF(N-i)(y(k))} }  • 

k 

Under  this  goal,  it  is,  for  example,  no  longer  efficient  for  the  manager  of  a  water 
supply  system  to  drain  the  reservoir  to  meet  current  demands  if  that  decision 
results  in  a  severe  shortage  in  the  next  stage.  To  solve  this  problem,  we  now 
need  an  estimate  of  F^. i>(y)  representing  the  expected  costs  of  one  following 
stage.  Assuming  that  optimal  decisions  will  also  be  applied  to  this  stage,  this 
expected  cost  function  is  simply  the  optimal  cost  function  of  the  previous  single 
period  solution;  i.e.. 


%-i)(yv-i)  =  J*(N)<?n) 

where  y#_i  =  .  If  we  again  determine  the  optimal  decisions  and  cost  for 

discrete  initial  states  xq  of  the  system,  we  can  specify  a  piece-wise  continuous 
optimum  control  function  u*(x)  and  associated  optimal  cost  function  J*(N. i)(x) 
for  a  two-stage  model. 

We  can  repeat  this  process  recursively  replacing  F{na)  by  F^. 2)  by  F^. 3) 
and  so  on,  until  we  incorporate  the  costs  of  all  N  stages  representing  the 
planning  horizon.  Each  recursion  steps  further  back  in  time  and  away  from  the 
zero-valued  expected  cost  function;  and  the  resulting  optimal  cost  function 
represents  expected  costs  for  a  planning  horizon  that  is  one  stage  longer. 

3.  Simple  Illustration 

The  method  is  clearly  illustrated  by  the  control  of  simple  water  supply 
system.  This  system  consists  of  a  single  reservoir  with  capacity  xmax  and  an 
uncorrelated  stream-flow  source.  The  state  of  the  system  is  described  by  the 
reservoir  level  x,  and  the  only  control  decision  u  is  how  much  to  release.  The 
transition  equation  is 
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y  =  x-u  +  w 


where  w  is  stream  flow  and  y  is  the  resulting  reservoir  level. 

Constraints  are  xmax  >  y  >  0  and  u>0.  In  standard  form  these  become 

u<x  +  w 

-U<Xmax-X-W 

and  -u<  0  . 

Assuming  that  planned  releases  can  only  use  previously  stored  water,  we  can 
ignore  stream  flow  in  the  first  constraint.  Also,  for  an  objective  considering  only 
water  supply,  the  second  constraint  can  be  dropped  since  spills  from  a  full 
reservoir  have  no  affect  on  optimal  control.  If  the  transition  function  produces  a 
final  reservoir  level  greater  than  the  maximum,  set  y  =  xmax . 

When  a  release  is  less  than  demand,  the  system  incurs  a  cost  C(u). 
Assuming  F^(y)  is  zero  for  any  ending  reservoir  level,  there  are  no  costs  beyond 
the  end  of  the  planning  horizon.  The  only  costs  that  affect  control  decisions  are 
those  accrued  between  the  current  and  last  stages. 

For  example,  if  there  is  only  one  remaining  stage,  the  efficient  release 
u*(x )  is  that  which  achieves 


=  min„{C(w)}  . 

Because  F(H)(y )  is  zero,  there  is  no  benefit  from  saving  water  in  the  reservoir  for 
future  use.  Thus  there  is  no  incentive  to  hedge.  The  efficient  decision,  u*,  is  to 
release  as  much  as  possible  to  meet  demand,  subject  to  the  constraint  u  <  x 
(Figure  II-6). 

If,  on  the  other  hand,  there  are  two  remaining  stages,  the  efficient  release 
becomes  that  which  achieves 

/V-i)(*)  =  mintt{C(M)  +  F(v.i)Cy)} 

subject  to  the  constraints  and  the  transition  equation.  Since  F^-i)  = 

=  min„-{C(w’)  +  mintt{C(u)}}  . 

There  now  is  an  incentive  to  hedge  in  the  current  period  if  the  cost  of  hedging  is 
less  than  the  benefit  of  the  extra  water  in  the  last  period.  As  expected  costs 
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incorporate  more  stages  and  a  longer  planning  horizon,  hedging  becomes  more 
common  (Figure  II-7). 

4.  Stochastic  Versus  Deterministic 

The  solution  of  the  above  problem  produces  an  expected  cost  function, 
F(y),  for  every  stage  of  a  planning  horizon.  We  can  think  of  these  expected  cost 
functions  as  the  sum  of  a  series  of  current  costs,  C,  resulting  from  a  control 
decision  applied  to  each  stage  remaining  until  the  end  of  the  planning  horizon.  If 
all  inputs  to  the  system  are  certain,  we  can  specify,  in  advance,  optimal  decisions 
and  states  of  the  system  by  a  single  control  schedule.  The  cost  function  is  a  sum 
of  known  costs 


N 

F(0)  =  X  {C(x,,u,)}  • 

t=  1 

However,  if  there  are  uncertain  inputs,  a  single  control  schedule  cannot  specify 
the  optimal  control  of  a  system  and  the  cost  function  is  a  sum  of  expected,  vice 
known,  costs. 

If  there  are  uncertain  inputs  w,  one  cannot  know  the  state  of  following 
stages  until  the  beginning  of  those  stages.  Thus  we  cannot  determine  optimal 
controls  and  costs  until  the  beginning  of  a  stage.  A  control  schedule  applied  to  a 
system  depends  on  unknown  future  states  of  the  system  resulting  from  realized 
values  of  w.  Because  the  state  at  any  particular  stage  is  unknown,  a  solution 
must  specify  optimal  controls  and  costs  for  any  possible  state. 

This  is  the  major  difficulty  in  solving  stochastic  dynamic  programming 
(SDP)  problems.  Because  we  must  specify  a  solution  for  every  possible  state  of  a 
system,  we  potentially  have  a  huge  number  of  calculations  to  perform  and  a 
huge  amount  of  data  to  store.  If  a  solution  can  be  specified  by  a  simple  function, 
the  parameters  of  the  function  specify  the  solution.  Generally,  however,  we 
cannot  specify  a  solution  by  a  simple  function.  Instead,  we  must  approximate 
the  solution  by  calculating  values  for  many  discrete  states — this  is  the  Discrete 
Dynamic  Programming  (DDP)  method.  To  maintain  a  particular  level  of  solution 
accuracy,  the  number  of  discrete  states  increases  exponentially  with  each 
additional  variable  added  to  the  state  vector.  Problems  with  more  than  two  or 
three  state  variables  generally  present  a  great  computational  challenge.  This  can 
be  addressed  partly  by  improved  optimization  methods,  such  as  discussed  in 
Chapter  IV.  Also  we  can  reduce  the  challenge  by  understanding  better  the 
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characteristics  of  a  solution,  thus  we  may  be  able  to  apply  simplified  models  and 
functions  that  reduce  the  number  of  discrete  states. 

Conjunctive-use  problems  necessarily  start  with  at  least  three  state 
variables,  and  this  is  when  the  model  that  is  generally  unrealistically  simple.  A 
conjunctive-use  model  requires,  at  a  minimum,  state  variables  to  represent  levels 
in  reservoir  and  aquifer  storage,  and  stream  flows.  Generally,  we  can  model  the 
state  of  a  reservoir  by  a  single  variable  for  storage,  though  each  reservoir  needs 
its  own  state  variable  unless  we  aggregate  their  characteristics.  On  the  other 
hand,  it  is  not  clear  that  the  state  of  aquifer  storage  could  ever  be  realistically 
modeled  by  a  single  state  variable  since  pumping  costs  depend  on  the  distributed 
characteristics  of  the  aquifer.  Stream  flows  are  usually  autocorrelated  for  the 
length  of  stages  practical  for  real-time  system  control,  and  state  variables  must  be 
added  to  incorporate  information  from  prior  flows  and  from  other  sources  such 
as  snow-pack  measurements.  For  example.  Table  1  outlines  the  state  variables 
that  may  be  required  to  model  the  EBMUD  system.  Because  of  the  number  of 
state  variables  required  to  model  realistic  conjunctive-use  problems,  it  is  unlikely 
that  such  problems  can  be  solved  by  traditional  SDP  methods. 

I  use  Gradient  Dynamic  Programming  (GDP)  of  Foufoula-Georgiou  and 
Kitanidis  [1988]  to  solve  the  problem.  GDP  is  a  new  method  that  allows  us  to 
solve  problems  with  a  greater  number  of  state  variables.  Chapter  IV  presents  the 
method.  The  algorithm  that  I  present  is  an  extension  of  the  original  method  to 
include  autocorrelation. 

5.  Convergence 

The  effects  of  current  decisions  are  preserved  and  transmitted  to  later 
stages  through  the  state  variables  y  that  result  from  the  decisions  and  inputs 
applied  to  the  initial  state  x.  F(y)  is  the  sum  of  all  expected  costs  starting  from 
state  y  and  applying  optimal  decisions  in  remaining  stages.  In  a  problem  with 
stochastic  inputs,  we  can  determine  F(y)  by  recursively  solving  problems  with 
increasing  numbers  of  stages  until  F( y)  represents  expected  costs  for  an  entire 
planning  horizon. 

The  cost  function  F(y)  influences  control  decisions  by  providing  the 
relative  costs  of  different  decisions.  The  absolute  value  of  F(y)  does  not 
influence  control  decisions,  thus  we  can  re-define  F(y)  as 


F(m)  =  /(<)*  -  C(t) 
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without  affecting  control  decisions.  Constant  q,)  standardizes  the  cost  function, 
allowing  us  to  more  easily  view  changes  with  time,  t  =  N,N-l,  N- 2, ...  1  ,  and 
avoiding  potentially  large-valued  functions.  Also,  with  an  appropriate 
definition,  we  can  use  C(t)  as  an  indication  of  the  overall  change  in  the  cost 
function,  such  as  the  extra  cost  expected  for  an  additional  stage  of  operation. 

Stedinger  et  al  [1984]  suggest  setting  c(0  to  the  average  gain  in  costs 
incurred  during  stage  t.  Another  nearly  equivalent  approach  is  to  set  cq)  to  the 
mean  or  other  measure  of  the  central  tendency  of  /(<)* .  A  third  approach  is  to  set 
C(j>  to  the  smallest  discrete  value  of  /(,)*  •  With  this  approach  F(m>  is  a  strictly 
positive  function;  though,  since  c(f)  is  based  on  a  single  discrete  value  instead  of 
an  average  of  all  discrete  value,  numerical  errors  can  be  more  significant. 

If,  after  N  periods,  the  assumed  cost  function  F(/v)( y)  has  a  negligible 
effect  on  current  decisions,  F(y),  u*(x),and  c  should  converge  (Figure  II-8). 

This  will  often  be  the  case  for  water  projects  where  the  planning  horizon  may  be 
tens  of  years.  In  these  cases,  we  can  solve  the  problem  with  fewer  recursions  by 
performing  only  as  many  as  required  to  obtain  convergence. 

When  convergence  is  achieved,  the  value  of  c  is  a  measure  of  stage  to 
stage  costs.  This  parameter  can  be  used  to  compare  performance  of  different 
system  configurations.  For  example,  this  parameter  can  be  used  to  compare  the 
benefit  of  adding  aquifer  storage  or  additional  reservoir  storage  to  a  system. 

Also,  the  current  state  of  the  system  has  only  a  transient  effect  on  the 
expected  cost  of  the  system;  if  we  look  far  enough  into  the  future,  we  expect  costs 
for  each  additional  period  of  operation  to  be  the  same  regardless  of  the  current 
state  of  the  system.  For  example,  though  the  system  may  currently  experience 
drought  or  flood,  the  expected  cost  of  future  operations  will  approach  the  long¬ 
term  average.  The  number  of  stages  required  to  achieve  convergence  indicates 
the  system’s  memory  of  control  decisions;  thus,  this  could  provide  some  insight 
into  differences  between  systems.  For  example,  limits  on  pumping  and  recharge 
rates  will  generally  result  in  a  slower  convergence.  This  is  related  to 
groundwater’s  relative  weakness  in  meeting  acute  short-term  water  supply  needs 
but  its  relative  strength  in  meeting  moderate  long  term  water  supply  needs. 

In  this  thesis,  c  is  the  smallest  discrete  value  of  /( .  I  use  c  only  as  a 
standardizing  constant  to  make  it  more  easy  to  view  changes  in  F(y);  we  prefer 
the  advantage  of  a  strictly  positive  F(y)  over  a  better  measure  of  change. 
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Recursive  Development  of  a  Control  Solution 


solution 

development 


time  - > 


.3  -2  -1  end 


max 


State  x 
of  the 
system 


min 


optimum  transitions- 


Fig.  II-5.  Real  time  of  decisions  progresses  to  the  right  while  development  of 
solution  recurses  to  the  left.  For  the  simple  water  supply  system  of  figure  1,  the 
sequence  of  optimum  transitions  are  those  that  might  be  expected  for  a  dry-wet- 
dry  sequence  of  stream  flows.  The  arrows  do  not  connect  because  we  develop  a 
cost-to-go  function  that  is  interpolated  between  the  discrete  values.  Thus  it  is  not 
necessary  to  require  that  transitions  end  on  the  discrete  states  used  to  span  the 
state  variables. 
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Optimal  Decisions  Applied  to  the  First  Recursion 


Optimal  decisions  minimize  current  and  future  costs: 

min  [fix) }  =  min{  C(n)  +  F(y ) }  . 

The  first  recursion  is  the  last  stage,  so  F(y )  is  zero  for  all  ending  states  y: 

min  {./(*) }  =  min{  C(u) } . 

For  each  discrete  reservoir  level  x we  determine  a  release  decision,  ut,  that 


This  results  in  an  expected  total  coster,)  and  change  in  reservoir  level: 


ftx)  =  C(u)  +  0, 
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Fig.  II-6.  Illustration  of  optimal  release  decisions  for  the  first  recursion  (last  stage) 
and  the  resulting  expected  costs  and  changes  in  reservoir  level. 
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Optimal  Decisions  Applied  to  the  Second  Recursion 


In  the  second  recursion,  the  expected  cost  function,  F(y),  is  no  longer  zero.  It  is 
now  equal  to  the  expected  total  cost  function  of  the  first  recursion.  Thus, 

min  { fix)  )  =  min  {  C(u)  +  F(y)  ) 

For  each  discrete  reservoir  level  jq,  we  again  determine  a  release  decision,  iq,  that 


This  again  results  in  an  expected  total  cost  fix  t)  and  change  in  reservoir  level,  but 
now  release  decisions  may  hedge  to  avoid  more  severe  rationing  in  the  following 


stage: 


fix)  =  C(m)  +  F(y), 
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Fig.  II-7.  Illustration  of  optimal  release  decisions  that  incorporate  hedging  when 
expected  future  costs  are  not  zero. 
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Number  of  Decision  Periods 
(prior  to  the  end) 


Fig.  II-8.  Increases  in  the  expected  costs  with  the  number  of  remaining  stages 
should  converge  to  a  constant  rate,  regardless  of  the  initial  state  of  the  system. 
Initial  states  xp  x2,  and  x3  illustrate  this,  where  x}  is  more  costly  than  x2  which  is 

more  costly  than  x3. 


E.  RESULTS 


The  solution  of  the  example  problem  demonstrates  four  important 
observations  about  the  optimal  control  of  conjunctive-use  systems:  (1)  efficient 
control  requires  effective  use  of  the  capabilities  of  both  reservoir  and  aquifer 
storage;  (2)  autocorrelation  significantly  affects  optimal  control;  (3)  analytic 
optimization  of  conjunctive-use  systems  is  important  and  perhaps  necessary  for 
achieving  optimal  control;  and  (4)  the  solution  of  the  conjunctive-use 
management  problem  has  characteristics  that  may  simplify  future  studies. 

1.  Efficient  Allocation 

The  solution  demonstrates  that  optimal  management  of  conjunctive-use 
systems  takes  advantage  of  the  strengths  of  both  mechanisms,  achieving  greater 
system  efficiency  and  reliability  than  alternate  systems  using  either  mechanism 
in  isolation.  Specifically:  (1)  optimal  management  results  in  lower  rationing 
costs  than  achievable  by  heuristic  management  (such  as  using  aquifer  storage 
only  as  a  back-up);  and  (2)  optimal  allocation  results  in  a  preference  for  aquifer 
storage  over  reservoir  storage  as  the  risk  of  near-term  rationing  declines. 

We  can  most  easily  observe  these  two  results  by  viewing  contours  of 
expected  rationing  cost.  Figure  II-9  displays  an  example  set  of  contours  plotted 
against  reservoir  and  aquifer  storage  levels.  These  contours  indicate  higher  costs 
for  lower  storage  levels.  Also,  the  surface  defined  by  these  contours  is  convex. 

Efficiency 

Given  a  total  amount  of  stored  water,  we  can  identify  a  minimum  cost 
allocation  between  reservoir  and  aquifer  storage.  A  general  principle  of 
economics  is  that  decisions  are  made  at  the  margin:  if  the  benefit  of  one 
additional  unit  of  groundwater  is  greater  than  the  benefit  of  one  additional  unit 
of  water  in  the  surface  reservoir,  we  would  prefer  to  recharge  groundwater  over 
filling  the  reservoir.  For  example,  the  allocation  of  300,000  acre-feet  of  water  is 
constrained  by  a  linear  function  A--A'  in  Figure  II-9.  Since  the  cost  curves  are 
convex,  we  can  easily  identify  a  minimum  cost  allocation.  For  300,000  acre-feet 
of  water,  the  minimum  cost  allocation  is  at  point  B. 

Given  expected  costs  as  in  Figure  II-9,  we  can  identify  the  lowest  cost 
allocation  using  both  reservoir  and  aquifer  storage.  However,  we  cannot  ensure 
that  this  will  be  achieved  by  control  decisions  because  of  uncertain  stream  flows. 
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The  precise  result  of  control  decisions  is  determined  by  values  of  stream  flow 
realized  in  the  future:  actual  flows  that  are  higher  than  expected  result  in  over 
allocation  to  the  reservoir,  actual  flows  that  are  lower  than  expected  result  in 
over  allocation  to  the  aquifer.  Also,  we  cannot  ensure  that  we  will  achieve  the 
lowest  cost  allocation  because  of  constraints  on  control  decisions  that  may  limit 
our  ability  to  reach  this  allocation. 

Though  we  cannot  ensure  that  the  lowest  cost  allocation  will  be  achieved 
by  control  decisions,  we  can  apply  control  decisions  that  attempt  to  reach  this 
allocation.  The  optimization  approach  presented  in  this  work  attempts  to  reach 
this  allocation,  subject  to  the  constraints  and  operational  costs  that  these 
decisions  entail.  In  contrast,  control  decisions  determined  heuristically  do  not 
recognize  the  optimal  allocation  and  thus  do  not  attempt  to  reach  it.  Thus, 
heuristic  control  rules  generally  are  less  efficient  than  rules  constructed  from  the 
solution  of  an  optimization  problem. 

Allocation 

Given  the  ability  to  identify  the  least  cost  allocation  of  water,  we  can 
determine  the  optimum  allocation  for  any  total  amount  of  stored  water.  Figure 
11-10  reorients  the  display  of  contours  in  Figure  II-9  to  display  more  clearly  the 
least  cost  allocation  of  any  total  amount  of  stored  water.  The  least  cost 
allocations  occur  at  the  trough  in  the  cost  contours  as  indicated  by  a  line  in 
Figure  11-10.  We  can  observe  the  least  cost  allocation  favors  storing  most  water 
as  groundwater,  except  when  storage  levels  are  very  low. 

These  contours  specify  costs  expected  in  November  when  stream  flow  in 
October  is  above  average.  The  least  cost  allocations  for  these  contours  favor 
groundwater  because  (1)  the  likelihood  of  near-term  rationing  is  low  in 
November  at  the  beginning  of  the  winter  rainy  season  and  (2)  the  likelihood  that 
future  stream  flows  will  spill  (and  be  lost)  from  a  full  reservoir  is  high  if  more 
water  is  stored  in  the  surface  reservoir.  Since  the  model  assumes  that  current 
demand  for  water  can  be  met  only  from  storage  and  not  from  current  stream 
flows,  it  is  necessary  to  store  a  small  amount  of  water  in  the  reservoir. 

2.  Effect  of  Autocorrelation 

When  stream  flow  time  series  are  autocorrelated,  a  below-average 
monthly  stream  flow  decreases  the  expected  value  of  later  stream  flows. 
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Expected  costs  increase  since  lower  expected  stream  flows  increase  the  likelihood 
of  rationing. 

Figures  11-11  displays  expected  cost  contours  for  November  when  stream 
flow  in  October  is  below  average.  Compared  to  Figure  11-10,  expected  costs  are 
greater  for  all  combinations  of  storage  levels  since  prior  stream  flow  is  20,000  vice 
60,000  acre-feet  for  October. 

Autocorrelation  influences  the  expected  costs  and  consequently  the 
optimal  decisions  that  allocate  water  to  users  and  between  storage  mechanism. 
Optimal  control  of  the  conjunctive-use  system  changes  decisions  that  allocate 
and  ration  water  for  users  and  that  allocate  water  between  reservoir  and  aquifer 
storage. 

Allocation  to  Users 

System  operators  must  continuously  update  control  decisions  based  on 
the  realized  and  expected  values  of  inputs.  With  each  decision,  system  operators 
make  control  decisions  that  balance  current  performance  against  uncertain  future 
performance.  Under  these  conditions,  rational  operators  will  hedge  their  control 
decisions  to  improve  future  performance. 

As  expected  costs  increase  with  decreasing  stream  flow,  rational  operators 
have  a  greater  incentive  to  hedge  by  rationing  the  water  supplied  to  users.  This 
rationing  reduces  the  possibility  of  more  severe  rationing  in  the  future. 

Figure  11-12  is  a  surface  plot  of  the  release  decision  for  various  storage 
levels.  The  figure  shows  that  less  water  is  released  as  levels  in  both  reservoirs 
decrease.  With  less  water  in  storage,  the  likelihood  of  future  shortages  increases, 
increasing  the  expected  costs.  Thus,  the  surface  plot  of  Figure  11-12  has  a  shape 
that  is  roughly  an  inverse  of  that  which  would  result  from  the  expected  costs  of 
Figure  11-11. 

Figure  11-12  displays  the  release  decisions  for  November  when  stream 
flow  in  October  is  20,000  acre-feet.  For  other  months  and  prior  stream  flow, 
releases  generally  decrease  with  increases  in  expected  cost,  and  vice-versa. 

Allocation  Between  Reservoir  and  Aquifer  Storage 

Besides  changing  decisions  that  ration  water  supplies,  stream-flow 
autocorrelation  also  influences  optimal  decisions  that  allocate  water  between  the 
dissimilar  mechanisms  of  reservoir  and  aquifer  storage.  Autocorrelation 
increases  the  time  that  prior  stream  flows  influence  future  flows;  however. 
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autocorrelation  also  decreases  the  severity  of  droughts  due  to  the  effective 
averaging  of  numerous  realizations  of  stream  flow.  Thus,  the  degree  of 
autocorrelation  in  a  stream-flow  time  series  significantly  affects  the  efficient 
allocation  of  stored  water  between  reservoir  and  aquifer  storage. 

As  we  saw  above.  Figures  11-10  and  11-11  display  the  expected  cost 
contours  for  November  when  prior  stream  flows  are  above  average  and  below 
average.  We  can  see  how  the  minimum  cost  allocations— indicated  by  lines 
connecting  the  minimum  cost  points  of  the  contours — change  with  prior  stream 
flow.  As  prior  stream  flows  decrease,  the  likelihood  of  near-term  water 
shortages  increases.  Comparing  Figure  11-11  to  Figure  11-10,  we  can  see  a 
reduced  preference  for  aquifer  storage  due  to  the  higher  risk  of  near-term 
rationing.  This  is  indicated  both  by  the  shift  in  the  minimum  cost  allocation  and 
by  the  flatter  cost  contours.  As  risk  of  near-term  shortage  increases,  it  is  less 
likely  that  the  reservoir  capacity  constraint  will  be  binding  on  decisions,  and  the 
minimum  expected  cost  is  not  sensitive  to  allocation.  When  the  risk  of  near-term 
shortage  is  high,  the  only  inefficient  allocations  are  those  that  make  exclusive  use 
of  aquifer  storage:  groundwater  can  meet  only  about  a  third  of  current  demand 

due  to  the  pumping  rate  constraint. 

The  minimum  cost  allocations  of  Figures  11-10  and  11-11  control  the 
groundwater  allocation  decisions  of  Figures  II-13a  and  II-13b.  Figure  II-13a 
displays  recharge  (positive  values)  and  pumping  (negative  values)  decisions  for 
November  when  stream  flow  in  October  is  above  average.  We  can  see  that 
optimal  control  prefers  groundwater  recharge  for  all  storage  levels  except  when 
the  reservoir  is  nearly  empty.  Comparing  Figure  II-13a  to  Figure  11-10,  we  can 
see  that  when  expected  costs  of  rationing  are  low,  we  should  store  water  in  the 
aquifer  to  meet  long-term  needs.  Conversely,  Figure  II-13b  displays  recharge 
and  pumping  decisions  for  November  when  stream  flow  in  October  is  below 
average.  We  can  see  that  optimal  control  prefers  groundwater  recharge  less  than 
in  Figure  II-13a.  The  differences  in  allocation  decisions  between  Figures  II-13a 
and  II-13b  reflect  the  differences  in  expected  costs  between  Figures  11-10  and  II- 
11. 

As  a  final  note.  Figures  II-13a  and  II-13b  display  a  dip  in  the  recharge 
decisions  as  groundwater  storage  becomes  large.  This  dip  may  be  a  result  of  the 
minor  pumping  and  recharge  costs  applied  to  the  objective  function.  The  dip 
may  also  be  an  artifact  of  the  optimization  method  or  an  error  in  the  computer 
code  used  to  obtain  the  solution:  I  have  observed  that  some  of  the  other  control 
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surfaces  have  anomalous  oscillations.  I  hope  that  subsequent  work  will  identify 
the  cause  of  these  oscillations. 

3.  Optimization  and  Simulation 

Burges  and  Maknoon  [1975]  note  that  "An  important  direction  for  research 
effort  would  be  to  examine  the  general  nature  of  'response  surfaces’  for 
conjunctive  use  problems.  If  they  turn  out  to  be  relatively  flat  then  it  may  be 
possible  to  use  the  simulation  approach  rather  than  direct  optimization 
procedures.  This  would  overcome  many  of  the  problems  of  dimensionality 
which  limit  the  use  of  optimization  models." 

As  we  can  see  from  the  above  results,  the  "response  surfaces"  for 
conjunctive  use  problems  are  sometimes  flat  and  sometimes  not.  When  the  risk 
of  water  shortages  is  high  in  the  example  problem,  the  expected  cost  contours  are 
flat  and  not  sensitive  to  allocation  between  reservoir  and  aquifer  storage. 
However,  the  example  problem  does  not  consider  the  costs  of  groundwater 
pumping  or  recharge.  Including  these  costs  will  tip  the  cost  contours  to  a  greater 
preference  for  reservoir  storage,  especially  when  the  risk  of  shortages  is  high. 

Because  of  the  different  characteristics  of  reservoir  and  aquifer  storage,  we 
can  conclude  that  the  "response  surfaces"  for  conjunctive-use  problems  are 
generally  not  "flat"— that  is,  the  values  of  the  objective  and  optimal  control 
decisions  change  with  the  state  of  the  system.  Thus,  optimization  of  conjunctive- 
use  problems  can  significantly  improve  the  control  of  conjunctive-use  systems. 
Under  these  conditions,  efficient  heuristic  control  of  conjunctive-use  systems  is 
unlikely.  However,  it  may  be  possible  to  specify  simple  functions  that  capture 
characteristics  of  the  response  surfaces,  thus  reducing  the  number  of  discrete 
states  that  must  be  considered. 

4.  Possible  Simplifications 

One  possible  simplification  that  could  be  applied  to  future  conjunctive-use 
problem  solutions  is  the  use  of  a  global  function  as  advocated  by  Buras  [1972]. 

We  can  observe  that  the  expected  costs  and  control  decisions  for  this  stochastic 
problem  are  quite  smooth.  This  is  largely  a  result  of  the  stochastic  inputs  which 
result  in  averaging  over  the  possible  input  realizations.  Because  of  this 
smoothness,  it  may  be  possible  to  describe  the  cost  and  control  functions  by  a 
global  polynomial  or  exponential  function.  Thus  we  could  avoid  some 
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discretization  of  the  state  space  and  solve  similar  problems  with  less 
computational  effort. 

Another  possible  simplification  is  the  use  of  more  constrained  operating 
policies  that  simplify  control  of  aquifer  storage.  Bogle  and  O’Sullivan  [1979] 
observe  that  optimal  operating  policies  derived  by  dynamic  programming  have  a 
"common  sense"  form  that  results  in  a  release  that  is  "minimum  for  storage  below 
some  level  and  maximum  for  storage  above  some  other  level."  Though  the 
results  presented  here  do  not  support  such  a  simplification,  there  may  be 
situations  in  which  an  operating  policy  presents  a  consistent  form  that  can  be 
used. 
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Expected  Cost  of  Water  Rationing 
Cost  Versus  Current  Storage 


Fig.  II-9.  Contour  plot  of  expected  costs  (in  arbitrary  units)  versus  amount  of 
stored  water  in  reservoir  and  in  aquifer. 

The  dashed  line  between  A  and  A'  indicates  possible  allocations  of  300,000  acre- 
feet  of  total  stored  water.  Point  A  represents  allocation  of  all  300,000  acre-feet  to 
reservoir  storage  with  an  expected  cost  of  520.  Similarly,  point  A'  represents 
allocation  of  all  300,000  acre-feet  to  aquifer  storage  with  an  expected  cost  of  460. 
Expected  cost  can  be  reduced  to  370  at  point  B  by  allocating  water  to  both 
reservoir  and  aquifer  storage. 
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Minimum  Cost  Allocation  of  Stored  Water 
Expected  Cost  Versus  Total  Amount  and  Allocation  of  Stored  Water 
(prior  month's  stream  flow  above  average) 


constant 

/  minimum 

— 800'"" 

cost 

^  cost 

contour 

allocation 

Fig.  11-10.  Contour  plot  of  expected  costs  versus  total  amount  and  allocation  of 
stored  water.  Amount  of  water  stored  in  reservoir  and  aquifer  can  be  measured 
diagonally  from  zero.  The  results  are  for  November  in  the  first  year  of  a  four-year 
planning  horizon  when  October's  stream  flow  was  above  average  (60,000  acre-feet). 

The  stippled  line  indicates  the  allocation  of  stored  water  that  produces  the  lowest 
expected  cost  from  rationing.  Point  B  marks  the  optimum  allocation  for  300,000 
acre-feet  determined  in  Figure  II-9. 
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Minimum  Cost  Allocation  of  Stored  Water 
Expected  Cost  Versus  Total  Amount  and  Allocation  of  Stored  Water 
(prior  month's  stream  flow  below  average) 
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Fig.  n-11.  Contour  plot  of  expected  costs  versus  total  amount  and  allocation  of 
stored  water.  The  results  are  for  November  in  the  first  year  of  a  four-year  planning 
horizon  when  October's  stream  flow  was  below  average  (20,000  acre-feet). 
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Supply  Decision 

Water  Supplied  Versus  Current  Reservoir  and  Aquifer  Storage  Levels 


Amount 

Supplied 


Reservoir 
Storage  Level 


Aquifer 
Storage  Level 


Fig.  11-12.  Surface  plot  of  water  supplied  to  meet  a  demand  of  58.3  per  month. 
Decisions  to  supply  less  than  this  amount  indicate  rationing.  Results  are  for 
November  in  the  first  year  of  a  four-year  planning  horizon  when  October's  stream 
flow  is  20.  Amounts  are  in  1000's  acre-feet. 
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Aquifer  Storage  Decision 
Amount  of  Water  Recharged  (+)  or  pumped  (-) 
Versus  Current  Reservoir  and  Aquifer  Storage  Levels 


Fig.  II-13a.  Surface  plot  of  change  in  aquifer  storage  level.  Positive  values 
indicate  recharge,  with  a  maximum  of  10  per  month.  Negative  values  indicate 
pumping,  with  a  minimum  of  20  per  month  (pumping  cannot  exceed  20,000 
acre-feet/month).  Results  are  for  November  in  the  first  year  of  a  four-year 
planning  horizon  when  October's  stream  flow  is  above  average  (60,000  acre-feet) 
Amounts  are  in  1000's  acre-feet. 


Fig.  II-13b.  Surface  plot  of  change  in  aquifer  storage  level.  Results  are  for 
November  when  October's  stream  flow  is  below  average  (20,000  acre-feet). 
Amounts  are  in  1000's  acre-feet. 


F.  CONCLUSIONS 


I  have  illustrated  the  optimal  control  of  a  conjunctive-use  problem  with 
uncertain  and  autocorrelated  inputs.  From  the  results,  we  can  observe  that 
differences  in  the  characteristics  of  reservoir  and  aquifer  storage  affect  optimal 
control.  This  is  especially  true  in  the  presence  of  inputs  that  are  autocorrelated. 

Results  from  an  analysis,  such  as  presented  here,  can  provide  practical 
information  to  planners  and  managers  of  water  supply  systems.  Planners  can 
anticipate  the  best  choice  of  mechanisms  to  enhance  the  capabilities  of  a  water 
supply  system  based  on  needs:  reservoir  storage  likely  will  be  a  better  choice  for 
short  term  or  seasonal  needs;  aquifer  storage  likely  will  be  better  for  inter-annual 
storage  needs.  Managers  can  obtain  efficient  operating  rules  that  are  applicable 
to  real  control  problems:  given  a  method  for  determining  the  optimal  control  of 
a  conjunctive-use  system,  control  decisions  can  be  reduced  to  simpler  rules 
useful  for  practical  operation  of  a  system. 

In  practice,  reservoirs  are  often  operated  according  to  rule  curves  that 
define  desired  storage  and  release  targets  [Young,  1967;  Loucks  et  ah,  1981]  that  are 
often  presented  in  graphical  form  as  guides  for  those  responsible  for  reservoir 
operation.  Many  rules  are  the  result  of  negotiation  between  competing  interests, 
and  may  be  legally  defined  to  protect  operators  and  managers  against  liability  or 
conflicts  among  beneficiaries  and  agencies.  Johnson  et  al.  [1991]  noted  that  most 
models  use  heuristic  guidelines  to  define  the  system's  operating  policy.  Often 
these  are  developed  for  application  to  a  specific  system;  however,  a  few  general 
rules  exist.  For  example,  the  "space  rule"  allocates  reservoir  storage  space  in 
proportion  to  expected  inflows  in  a  parallel  water  supply  system.  The  similar 
"New  York  City  rule"  allocates  reservoir  storage  space  to  obtain  an  equal  spill 
probability  for  each  reservoir  [Johnson  et  al,  1991]. 

Though  optimization  efforts  can  improve  our  management  of  conjunctive- 
use  systems,  they  are  still  of  limited  practical  benefit.  Economic,  social,  and  legal 
aspects  of  the  conjunctive-use  problem  are  still  absent;  as  Burges  and  Maknoon 
[1975]  state  "There  is  no  advantage  in  building  or  using  a  model  of  a  conjunctive 
ground-surface  water  system  that  includes  considerable  hydrologic  detail  but 
neglects  to  adequately  represent  legal  and  economic  nuances." 

Also,  Yeh  [1985]  notes  that  few  practical  rules  have  yet  been  developed 
through  application  of  optimization  techniques  developed  in  the  literature: 
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"reservoir  operators  are  still  very  reluctant  to  use  optimization  models  for  their 
day  to  day  scheduling  of  water  releases  and  power  generations."  Often  this  is 
because  (1)  operators  have  not  been  directly  involved  in  the  development  of  the 
computer  model,  (2)  most  published  papers  deal  with  simplified  reservoir 
systems  and  are  difficult  to  adapt;  and  (3)  there  are  institutional  constraints  that 
impede  user-research  interactions. 

In  addition,  though  conjunctive  use  potentially  can  improve  the  efficiency 
and  reliability  of  a  water  supply  system,  stumbling  blocks  remain  in  its 
implementation.  Lettenmaier  and  Burges  [1979]  point  out  that  practical 
implementation  of  aquifer  storage  is  inhibited  by  the  manner  in  which  costs  of 
surface  storage  facilities  are  allocated,  and  that  legal  questions  about  the  control 
of  subsurface  storage  remain  to  be  answered. 
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III.  Modeling  Issues 


In  any  modeling  effort,  we  generally  must  make  simplifying  assumptions 
about  system  characteristics.  These  are  required  to  develop  a  problem  that  we 
can  solve.  For  dynamic  systems,  the  modeling  of  uncertainty  and  the  stage-wise 
evolution  of  the  system  present  particular  challenges.  In  addition,  we  also  must 
consider  factors  that  affect  the  efficiency  and  tractability  of  the  solution  algorithm 
used  to  solve  a  conjunctive-use  problem. 

A.  QUANTIFYING  SYSTEM  CHARACTERISTICS 

The  general  application  of  optimization  techniques  to  water  management 
problems  has  resulted  in  some  topics  of  controversy.  Here  I  consider  a  few  of 
these  topics  and  how  they  should  influence  our  optimization  efforts. 

In  practice,  there  are  a  number  of  potential  difficulties  that  we  must 
consider  when  modeling  water  systems.  Two  difficulties  have  been  especially 
controversial:  (1)  quantifying  the  costs  associated  with  management  decisions; 
and  (2)  characterizing  the  probability  distributions  of  stochastic  components  in 
the  system.  In  addition  to  these  two  difficulties,  we  also  must  select  state 
variables  cautiously  to  ensure  a  tractable  problem  that  is  still  sufficiently  realistic 
to  be  of  practical  benefit. 

1.  The  Objective  Function 

Often  managers  must  consider  costs  which  are  poorly  defined.  Estimated 
costs  of  rationing  and  floods  are  approximate  at  best,  and  may  neglect  important 
impacts.  Estimated  costs  for  other  values  may  be  entirely  absent:  costs 
attributed  to  environmental  effects,  other  third-party  effects  (i.e.,  externalities), 
risk  aversion,  recreation,  aesthetics,  etc. 

Nevertheless,  though  these  factors  are  difficult  to  quantify,  there  generally 
are  clear  qualitative  differences.  Though  we  may  not  know  the  precise  costs 
associated  severe  flooding,  we  know  that  costs  increase  with  the  severity. 

To  obtain  quantitative  measures  of  cost,  we  often  use  values  provided  by 
market  transactions.  For  example,  revenue  obtained  from  the  sale  of 
hydropower  may  be  an  appropriate  measure  of  benefit.  However,  though 
market  values  may  be  convenient  and  widely  accepted,  they  may  be 
inappropriate  when  markets  do  not  incorporate  important  externalities. 
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No  Specified  Costs 

Yao  and  Georgakakos,  [1993]  observe  that  "Reservoir  systems  are  usually 
multiobjective  and  their  operation  involves  a  wide  spectrum  of  considerations 
including  interagency  coordination,  legal  and  regulatory  issues,  and  public 
pressures,  among  others.  Such  criteria  cannot  be  modeled  and,  therefore,  an 
operation  which  maximizes  or  minimizes  one  particular  objective  is  meaningless 
and  is  not  practiced."  Instead,  they  advocate  a  set  approach  that  successively 
limits  available  decisions  in  order  to  maintain  the  system  within  constraints. 

However,  Yao  and  Georgakakos  recognize  that  their  approach  is  limited, 
primarily  because  it  fails  to  distinguish  the  desirability  of  different  feasible 
solutions.  Also,  for  a  long  enough  time  horizon,  they  state  that  there  may  be  no 
feasible  solution.  While  their  approach  can  define  the  effect  that  different 
constraints  have  on  the  remaining  management  options  (and  therefore  can  be 
useful  to  policy  makers  in  determining  the  constraints  to  apply  to  operation  of  a 
system),  it  cannot  tell  a  manager  what  the  optimum  control  strategy  should  be. 
Nor  can  it  tell  the  operator  what  real-time  control  decisions  to  implement.  In 
addition,  constraints  used  by  the  set  approach  to  establish  "feasible"  controls  are 
not  truly  bounds  on  feasible  solutions;  often  they  are  merely  a  result  of 
management  decisions  or  of  negotiation. 

Hybrid  Specification  of  Costs 

In  contrast  to  the  work  of  Yao  and  Georgakakos,  most  efforts  to  optimize 
water  system  management  assume  that  an  objective  function  can  be  developed  to 
determine  optimum  control  strategies.  Often,  however,  these  efforts  use  a  hybrid 
approach  that  avoids  a  complete  definition  of  costs. 

When  confronted  with  poorly  defined  costs  or  other  difficult  to  model 
criteria,  it  has  been  common  for  model  developers  to  use  the  extreme  assumption 
that  nothing  is  known  about  these  costs  or  effects.  Often  this  approach  is  taken 
to  avoid  an  explicit  statement  of  costs,  such  as  that  associated  with  the  loss  of 
human  life,  because  of  the  ethical  and  political  problems  associated  with  such 
statements. 

However,  management  decisions  always  contain  an  implicit  definition  of 
costs,  even  for  those  criteria  for  which  we  have  not  directly  specified  costs.  For 
example,  chance  constraints  are  often  used  when  there  is  significant  uncertainty 
in  the  costs  or  probabilities  associated  with  extreme  events.  In  flood  control 
problems,  the  use  of  chance  constraints  allows  us  to  avoid  quantifying  the  costs 
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associated  flood  events  that  exceed  some  threshold  of  chance.  This,  however, 
ignores  our  knowledge  that  larger  floods  cause  greater  damage,  in  spite  of 
exceeding  a  threshold  of  chance.  The  use  of  chance  constraints  results  in  an 
inferior  objective  value  that  considers  the  cost  of  all  extreme  events  to  be  the 
same.  Similarly,  the  solution  boundaries  established  by  the  set  approach  of  Yao 
and  Georgakakos,  [1993]  are  equivalent  to  specifying  an  objective  which  has  a  zero 
cost  for  "feasible"  solutions  and  one  which  has  very  large  costs  for  "infeasible" 
solutions. 

Specified  Costs 

Any  approach  that  ignores  information  or  incorporates  it  in  the  objective 
by  a  prescribed  approach  will  generally  perform  more  poorly  than  one  that 
makes  appropriate  use  of  this  information,  no  matter  how  poorly  quantified. 

On  the  other  hand,  some  methods  have  been  developed  to  provide  a 
measure  of  non-market  costs  and  benefits.  One  of  the  more  popular  is  the 
contingent  valuation  method  (for  discussion  and  a  recent  application,  see 
Whittington  et  al,  [1993]).  Though  these  other  methods  can  contribute  to  the 
solution  of  the  conjunctive-use  problem  presented  in  this  thesis,  I  assume  that  all 
costs  associated  with  system  operations  are  quantifiable  and  that  the  cost 
function  used  to  evaluate  the  management  solution  is  known.  The  difficulties 
associated  with  including  difficult  to  quantify  (often  non-market)  costs  is  not 
considered  here. 

2.  The  Probability  Distribution  of  Stochastic  Inputs 

Describing  the  continuous  probability  distributions  of  stochastic 
components,  such  as  stream  flow,  is  the  second  difficulty  encountered  in 
modeling  efforts. 

In  many  systems,  the  amount  of  information  is  insufficient  to  determine  a 
unique  model  for  a  probability  distribution.  The  best  model  will  often  depend 
on  the  problem  under  consideration.  For  certain  types  of  problems,  the  reliability 
of  a  control  solution  may  be  poor.  Even  for  the  same  data  set,  such  as  historical 
stream  flows  of  a  particular  river,  a  control  solution  for  extreme  events  such  as 
floods  will  be  less  reliable  than  one  for  mean  characteristics  such  as  water 
supply.  The  probability  of  extreme  events,  such  as  of  floods  and  droughts,  is 
almost  never  known  with  much  accuracy.  Even  our  knowledge  of  the  mean  and 
variance  of  many  stochastic  elements  can  be  very  poor. 
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Because  of  this  difficulty,  some  authors  have  advocated  approaches  to 
uncertainty  that  may  use  less  knowledge  of  a  distribution  than  we  have  available 
[Yao  and  Georgakakos,  1993].  The  problem  with  these  approaches  is  similar  to  that 
of  quantifying  an  objective  function;  using  less  knowledge  does  not  improve  our 
solutions,  it  merely  avoids  decisions  that  may  be  awkward.  Others  have 
suggested  methods  to  creatively  incorporate  available  information  without 
developing  a  continuous  model  of  a  distribution  [Kelman  et  al.,  1990].  Even  if  we 
decide  that  it  is  best  to  use  a  probability  model,  selecting  an  appropriate  model 
can  be  difficult  [Vogel  and  Fennessey,  1993]. 

For  problems  in  which  managers  are  still  unwilling  or  unable  to  quantify 
certain  costs  or  distributions,  various  approaches  exist.  Foufoula-Georgiou  and 
Kitanidis  [1988]  use  chance  constraints  with  their  GDP  method.  Kelman  et  al. 

[1990]  avoids  using  a  model  of  flow  distribution  by  a  sampling  stochastic 
dynamic  programming  method  that  preserves  many  of  the  important 
characteristics  of  a  distribution.  Yao  and  Georgakakos,  [1993]  use  a  set  approach 
that  avoids  the  use  of  both  an  explicit  objective  and  a  model  of  flow  distribution, 
though  it  cannot  be  used  for  real-time  control  except  perhaps  during  emergency 
operations. 

3.  Incorporating  Uncertainty 

There  are  many  sources  of  uncertainty  in  system  modeling.  I  discuss  here 
two  sources  of  uncertainty  that  can  be  modeled  by  using  stochastic  variables. 
Both  affect  the  transition  of  the  system  from  one  period  to  the  next.  The  first  is 
uncertainty  in  system  inputs,  such  as  stream  flow.  The  second  is  uncertainty  in 
the  values  of  state  variables  used  to  describe  the  system.  Though  both  are  a 
result  of  limited  knowledge  about  the  system,  the  first  is  generally  attributed  to 
the  stochastic,  uncontrollable  nature  of  input  process,  and  the  second  is  generally 
attributed  to  uncertainty  in  measurements  or  decisions  over  which  some  control 
exists. 

Three  cases  of  uncertainty  are  discussed  below  followed  by  an  illustration 
of  how  they  might  be  included  in  a  system  model.  The  first  two  deal  with 
stochastic  uncontrollable  inputs,  one  that  results  from  a  random  uncorrelated 
time  series  and  one  from  an  autocorrelated  series.  The  third  deals  with  an 
uncertain  state  variable  where  uncertainty  can  be  reduced  by  more  precise 
measurements. 
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Based  on  the  previously  discussed  mathematical  form  of  the  model, 
uncertain  inputs  and  uncertain  measures  of  a  system's  characteristics  affect 
control  decisions  by  making  the  end-of-stage  value  of  state  variables  uncertain. 
Optimum  decisions  minimize  the  sum  of  current  and  future  costs,  and  future 
costs  depend  on  the  values  of  y.  Again,  the  mathematical  form  of  the  control 
problem  is 

F(t)(x)  =  minu  {  C(x,u)  +  £  ^/(r+iM*))}  1 

k 

where  y^  =  T(x,u,w^)  is  a  function  of  the  current  state  x,  decisions  u,  and 
stochastic  variables  w. 

In  the  first  case,  the  state  of  the  system  is  influenced  by  a  random 
uncorrelated  input,  a  stochastic  variable  w-  is  added  to  the  vector  w  to 

represent  the  realization  of  this  input  during  the  current  stage.  Though  the  value 
of  Wj  is  uncertain,  its  probability  distribution  is  assumed  known  and  capable  of 

being  modeled.  Since  the  input  is  uncorrelated  with  prior  realizations  or  other 
characteristics  of  the  system,  no  additional  state  variables  need  be  added  to  x. 

In  the  second  case,  the  state  of  the  system  is  influenced  by  a  random  input 
that  is  autocorrelated,  one  or  more  state  variables  must  also  be  added.  This  is 
necessary  because  the  realized  values  of  past  inputs  influence  optimal  decisions. 
The  added  state  variables  represent  prior  realizations  of  the  input;  however,  if 
the  model  of  autocorrelation  requires  many  prior  stream-flow  values,  much  of 
the  information  could  be  incorporated  by  a  single  state  variable  representing  a 
single  stream-flow  forecast  [Stedinger  et  al,  1984].  These  state  variables  condition 
the  distribution  of  the  random  variable,  and  discrete  realizations  used  to  span 
this  distribution  are  specified  as  functions  of  these  state  variables.  For  certain 
simple  models  of  probability  distribution  and  autocorrelation,  the  effect  of  the 
state  variables  on  the  stochastic  variable  could  be  built  into  the  transfer  equation 
and  the  stochastic  variable  used  to  model  only  the  independent  residual 
variation. 

In  the  third  case,  the  value  of  a  state  variable  has  a  significant  amount  of 
uncertainty,  and  a  stochastic  variable  can  be  included  to  incorporate  this 
uncertainty.  Though  not  a  physical  input  to  the  system,  the  stochastic  variable 
represents  additional  uncertainty  in  identifying  the  effect  of  control  decisions  on 
the  state  of  a  system.  If  this  uncertainty  is  due  to  measurement  error,  we  may 
have  some  control  over  this  uncertainty  by  making  additional  measurements.  In 
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a  dynamic  system,  measurements  made  over  time  present  an  uncertain  time 
series  that  may  also  be  autocorrelated  or  may  otherwise  depend  on  state 
variables.  For  example,  if  measurement  errors  of  a  state  variable  change  with  the 
variable's  magnitude,  then  the  stochastic  variable  representing  this  error  depends 
on  the  state  variable. 

In  each  of  the  above  cases,  adding  stochastic  variables  allows  the 
specification  of  a  transfer  function  that  has  a  precise  mathematical  form  similar  to 
that  discussed  in  Chapter  II.  Assuming  the  system  is  linear,  the  transfer  equation 
can  be  modeled  by  a  system  of  equations 


y  =  T(x,u,w)  =  Ax  +  Bu  +  Cw 


[Yao  and  Georgakakos,  1993].  Suppose  a  water  supply  system  is  described  by  two 
state  variables  where  xl  =  known  reservoir  level,  x2  =  either  the  prior  month's 
stream  flow  or  a  snowpack  measurement  ( Kelman  et  al,  1990).  If  x2  is  the  prior 
month’s  stream  flow,  we  can  incorporate  autocorrelation  by  a  first-order  auto¬ 
regressive  model.  If  x2  is  a  measurement  of  snowpack,  we  can  obtain  stream- 
flow  estimates  based  on  the  correlation  between  snowpack  and  snow-melt 
runoff.  Either  of  these  values  is  subject  to  measurement  error,  but  snowpack 
measurements  are  much  more  imprecise  and  subject  to  revision  by  later 
measurements.  Given  that  there  is  one  control  decision  for  this  system,  iq  = 
water  release  to  meet  current  demand,  the  resulting  state  y  depends  on  iq  and 
Wj  =  WjOq)  =  current  month's  uncertain  stream  flow. 

Consider  first  the  model  with  x2  =  the  prior  month's  stream  flow.  The 
transfer  equation  can  be  modeled  as 
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In  this  case,  the  uncertain  input  vv1  affects  both  end-of-stage  state  variables,  y2 
and  y2.  Also,  the  dependence  of  Wj  on  x2  could  be  model  as  first-order  auto¬ 
regressive,  where 

vq  =  jj.  +  pix2-/j)  +  eoV  1  -p2  . 

Because  of  the  simple  model  of  stream  flow  autocorrelation,  the  dependence  of 
current  stream  flow  on  prior  realizations  could  be  included  in  the  transfer 
function  as 
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y  =  T(x,u,w  )  = 
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where  p  is  the  sample  lag-one  autocorrelation  coefficient.  Using  this  model,  wx 
is  now  independent  of  x.  A  problem  modeled  with  this  transfer  function  could 
be  solved  by  the  algorithm  of  Foufoula-Georgiou  and  Kitanidis  [1988]  without 
extension. 

Alternately,  consider  the  model  with  x2  =  a  measurement  of  snowpack. 
Vector  w  now  contains  two  stochastic  variables:  Wj  =  current  month's 
uncertain  stream  flow,  and  w2  =  random  error  in  the  season's  stream-flow 
estimate.  For  this  system  the  transfer  equation  can  be  modeled  as 
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In  this  case,  the  uncertain  input  Wj  still  affects  both  end-of-stage  state  variables, 
y1  and  y2.  Also,  Wj  still  depends  on  the  state  of  the  system  as  some  function  of 
x2.  In  addition,  the  random  error  w2  now  affects  the  value  of  state  variable  y2, 
but  is  independent  of  the  values  of  the  state  variables  x. 

Both  the  prior  realization  of  stream  flow  and  the  stream-flow  estimate  for 
the  snow-melt  season  could  be  used  in  the  prediction  of  current  stream  flow; 
however,  this  would  be  at  the  cost  of  an  added  dimension  to  the  state  vector.  For 
example,  let  x1  =  known  reservoir  level,  x2  =  the  previous  month's  stream  flow, 
and  x3  =  a  measurement  of  snowpack.  Vector  w  is  defined  as  in  the  previous 
paragraph.  The  transfer  equation  can  be  modeled  as 
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In  this  case,  the  distribution  of  Wj  now  depends  on  a  model  that  includes  both 
x2  and  x3. 


4.  Modeling  Stochastic  Variables 

The  distribution  of  each  stochastic  variable  is  modeled  by  choosing 
discrete  realizations  w(yt)  of  the  multi- variate  distribution  w,  and  with  each  w(jk) 
we  can  specify  an  associated  probability.  For  example,  if  a  stochastic  variable  Wj 
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is  normally  distributed,  its  distribution  could  be  approximated  by  discrete  values 

wm  at 

G  =  0,  ±1,±2, . . . 

standard  deviations  from  the  mean.  These  realizations  could  then  be  used  to 
approximate  the  distribution  over  the  intervals 

....  [-2.5, -1.5],  [-1.5, -0.5],  [-0.5, 0.5],  [0.5, 1.5],  [1. 5,2.5], . . . 


with  probabilities 


Pn  —  0.383,  P[i± i g  —  0.242.  Pji+on  —  0.061,  ... 


which  are  the  weights  applied  to  each  realization  when  calculating  expected 
values.  Derivatives  dw^dx  are  estimated  by  changes  in  the  k' th  ordered 
discrete  realization  while  maintaining  the  same  deviation  from  the  mean. 
Given  some  simplifying  assumptions  discussed  below,  this  matrix  could  also  be 
approximated  by  the  covariance  matrix  Cov[w,x]  . 


Determining  Derivatives 

It  may  be  possible  to  calculate  the  elements  of  this  matrix  analytically 
using  the  probability  distribution  model  for  w.  In  the  model  developed  in 
Chapter  II,  I  was  able  to  use  the  formula 


c|w 

9x 


0, 


0, 


wl-*r 


based  on  the  three-parameter  log-normal  model  developed. 

In  some  situations  analytic  determination  of  derivatives  may  be  difficult. 
As  an  alternative,  each  value  dwj  ^dxj.  ^  can  be  approximated  by  regressing  Wj 
on  Xj.  from  historical  data.  Assuming  that  changes  in  Xj-  have  the  effect  of 
shifting  the  distribution  of  Wj  without  changing  the  shape  (an  appropriate 

simplifying  assumption  with  the  amount  of  data  normally  available), 
dwj^ftXj-  ^  will  be  equal  to  the  slope  of  the  regression  equation  at  x^ 
regardless  of  realization  (ik)  (Figure  III-l).  If  historical  data  are  limited,  only  a 
simple  linear  regression  of  Wj  on  x}-  may  be  possible.  In  this  situation, 
dwj  (kftXj"  (t)  will  be  constant  for  any  x^  ^  and  any  Wj^ .  Matrix  dw/dx  is  then 
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a  constant  matrix  R  which  could  be  calculated  just  once  at  the  beginning  of  the 
problem. 

Such  a  matrix  R  could  be  calculated  from  covariances  as  follows  [Draper 
and  Smith,  1981,  Chapter  2] 

rjf  =  Cov(wjjf  )VVar(wy)/VVar(x/) 

and  is  related  to  the  matrix  P  of  correlation  coefficients  between  w  and  x 
where  the  terms  of  P  are  calculated  as  follows 

pjj-  =  Co viwjjCj,  )/JVai(Wj)  Var(.Xy.)  . 

Thus,  the  matrix  8w/5x  preserves  the  correlation  between  w  and  x  and  the 
effect  of  predictors  on  the  stochastic  variables. 

Other  alternatives  to  the  above  approach  exist.  If  correlation  between 
stochastic  inputs  and  state  variables  can  be  explicitly  modeled,  an  alternate 
approach  may  be  to  incorporate  this  knowledge  directly  in  the  transition 
function.  This  leaves  only  uncorrelated  residual  uncertainty  to  be  modeled  by 
stochastic  variables  w,and  dw/dx  =0. 

5.  Incorporating  Trends,  Discount  Rates,  and  Seasonality 

It  is  possible  to  also  incorporate  various  other  system  characteristics  into 
this  approach,  such  as  trends,  discount  rates,  and  seasonality.  These 
characteristics  could  produce  changes  in  the  cost  function,  transition  function,  or 
constraints. 

Trends  can  be  incorporated  by  allowing  the  current  cost  function  to 
change  with  time,  though  the  change  would  have  to  be  applied  in  the  reverse  of 
its  actual  development  since  the  recursions  go  backwards  in  time.  The  recursive 
update  of  the  expected  cost  function  would  then  take  the  form 

F(Jx)  =  minu  {  C(t)(x,u)  +  X  lPf(t+i  )<?(*)»  ) 

k 

where  C,t)  varies  with  each  recursion.  Similarly  a  trend  could  be  applied 
through  changes  in  a  constraint  applied  to  each  recursive  solution.  Such  an 
approach  would  allow,  for  example,  changing  water  availability  resulting  from 
use  of  senior  reserved  water  rights,  such  as  water  reserved  for  area-of-origin  or 
for  federal  lands. 
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A  discount  rate  of  a  percent  can  be  applied  by  multiplying  the  updated 
expected  cost  function  by  the  value  (5  =  100  /  (100  -  a)  in  each  recursion  to 
reduce  the  present  value  of  future  costs: 

F(0(x)  =  minu  {  C(x,u)  +  0  £  [pkF(t+  ip^)}  }  . 

k 

The  effect  of  a  discount  rate  applied  to  the  entire  planning  period  can  be  more 
easily  observed  in  a  deterministic  problem  where  the  cumulative  effect  of  a 
discount  rate  over  the  entire  planning  horizon  would  be  represented  by 

N  ,  N  t 

F(0)(y)  =  minu  (  X  P  [C(X(, ),“(,))]  }  =  X  p  [C(X(, ),“(,)*)]  . 

<= 1  /=1 

In  a  deterministic  problem,  [u(1)*, . . .  u^*]  is  the  single  control  schedule 

needed  to  specify  the  solution  for  all  stages. 

Seasonality  can  be  incorporated  by  a  variety  of  approaches,  such  as  by 
developing  a  unique  expected  cost  function  for  each  season  or  by  adding 
another  variable  to  the  state  vector.  While  the  use  of  a  state  variable  is  more 
consistent  with  the  definition  of  the  state  of  a  system,  it  requires  development  of 
a  variable  whose  value  is  periodic.  Oscillating  functions  such  as  a  sine  or  cosine 
could  be  used,  but  they  may  imply  intra-seasonal  relationships  that  do  not  exist. 
Seasonality  is  most  often  incorporated  by  developing  a  unique  expected  cost 
function  for  each  season.  The  backwards  stage-by-stage  solution  of  an  expected 
cost  function  develops  as  previously  discussed,  but  costs,  constraints,  stochastic 
distributions,  or  even  the  transition  function  may  vary  periodically  with  the 
current  stage's  season.  We  can  represent  this  by 

F(t)(x)  =  minu  {C(T)(x,u)+  £  {Pk,r+iF(t+i)iy(k))}  } 

k 

subject  to  constraints  of  form  g(l^u  <  d(T)(x) 

where  t  is  an  index  for  the  season  and  is  specified  uniquely  for  any  stage  t. 
Note  that  F(I)  =  minu  {/(/+i) }  regardless  of  the  season.  For  sufficiently  long 
planning  horizons,  control  solutions  still  converge,  but  only  when  solutions  for 
equivalent  seasons  are  compared.  For  example,  if  a  monthly  decision  period  is 
used  as  in  Chapter  II,  convergence  of  control  decisions  is  observed  when 
comparing  the  most  recent  twelve  recursions  with  the  previous  twelve. 
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Regression  of  Stochastic  Variables  on  State  Variables 


6: 


B.  COMPROMISES 


In  addition  to  difficulties  in  quantifying  system  characteristics,  modeling 
also  requires  compromises  between  accuracy  and  feasibility.  Here  I  present  two 
characteristics  of  the  conjunctive-use  problem  that  required  consideration  in 
developing  the  problem  presented  in  Chapter  II. 

1.  Selecting  State  Variables 

As  Stedinger  et  al.  [1984]  point  out,  the  better  the  hydrologic  state  variables 
selected  to  predict  stochastic  characteristics  of  the  system,  the  better  control 
decisions  will  be.  However,  as  they  also  point  out,  if  all  information  on  the  state 
of  the  system  is  included,  the  dimension  of  x  can  become  large.  It  is  necessary 
for  us  to  choose  predictors  that  transmit  as  much  information  about  a  stochastic 
process  as  possible  with  the  fewest  number  of  state  variables  possible. 

This  can  often  be  accomplished  by  representing  several  characteristics  of 
the  system  by  a  single  variable.  For  example,  Stedinger  et  al.  [1984]  suggest  that 
an  estimate  of  future  stream  flow  can  be  useful  in  summarizing  information 
about  prior  flows.  Kelman  et  al.  [1990]  use  an  estimate  of  the  snow-melt  season's 
runoff  as  a  state  variable,  switching  to  the  prior  month's  flow  when  this  provides 
better  information  (and  even  dropping  the  extra  state  variable  entirely  for  those 
months  when  month-to-month  flow  correlations  were  modest). 

Incorporating  Autocorrelation 

When  prediction  is  based  primarily  on  autocorrelation,  traditional  auto¬ 
regressive  and  moving  average  models  can  be  employed.  A  practical  description 
of  observed  autocorrelation  will  rarely  require  the  addition  of  more  than  one  or 
two  state  variables.  Chat  field  [1989]  provides  a  good  practical  guide  in  fitting  a 
model,  and  Box  and  Jenkins  [1976]  provide  an  in  depth  analysis  of  the  popular 
ARIMA  model. 

Another  approach  to  including  the  information  of  predictors  without 
using  extra  state  variables  is  the  adaptive-control  method  [Bras  et  al,  1983; 
Stedinger  et  al,  1984].  This  approach  does  not  include  all  predictors  in  the  state 
vector  but  instead  accounts  for  the  effect  of  these  predictors  by  making  transient, 
non-stationary  modifications  of  distributions  of  the  stochastic  variables.  A  major 
advantage  of  using  adaptive  control  is  that  it  reduces  the  number  of  state 
variables  needed  to  describe  the  state  of  the  system.  A  main  disadvantage, 
however,  is  that  it  does  not  allow  development  of  a  global  solution.  Adaptive- 
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control  develops  only  the  solution  which  is  appropriate  for  a  single  state  of  the 
system  and  for  particular  values  of  the  predictors.  Thus  a  new  solution  must  be 
developed  for  each  period  and  each  change  in  the  predictors. 

2.  Length  of  a  Stage 

Stages  are  increments  in  time  or  space  identified  as  a  unit  for  the 
application  of  control  decisions.  For  example,  a  real-time  management  problem 
requires  that  control  decisions  be  updated  periodically;  a  process  control 
problem  requires  that  decisions  be  made  for  each  section  of  the  assembly  line  or 
treatment  path.  Generally,  a  stage  is  defined  as  the  maximum  time  that  a  system 
can  be  operated  without  updating  control  decisions.  Stages  can  also  identify 
spatial  lengths  such  as  when  optimizing  the  steady-state  degradation  of  effluent 
down  a  stream  channel  [Cardwell  and  Ellis,  1993]. 

Each  stage  is  represented  by  one  recursion  of  the  solution.  The 
appropriate  length  of  a  stage  varies  with  the  problem  being  considered. 
Generally,  short  stages  are  desirable  for  more  accuracy;  however,  there  are  often 
practical  limits  on  how  short  a  stage  can  be,  and  the  length  of  a  stage  is  a 
compromise  between  accuracy  and  practicality. 

The  Improved  Efficiency  of  Short  Stages 

The  length  of  a  stage  affects  solution  accuracy  in  several  ways.  Stochastic 
inputs  to  the  system  make  the  cost  of  decisions  uncertain.  Inaccuracies 
introduced  by  this  uncertainty  can  frequently  be  reduced  by  choosing  a  short 
stage.  As  shorter  stages  are  chosen,  current  costs  C(x,u)  become  negligibly 
different  from  the  expected  value  Ew{C(x,u,w)}  .  This  is  when  inputs  are  highly 
autocorrelated. 

For  example,  in  a  water  supply  system  with  a  random  demand 
component,  the  actual  loss  incurred  from  a  mismatch  between  demand  and  a 
release  decision  will  not  be  precisely  known  until  the  end  of  the  stage.  However, 
because  demand  for  water  is  highly  autocorrelated — especially  on  short  time 
scales— the  shorter  the  stage  used,  the  smaller  the  error  in  the  cost  estimate. 

Another  related  consideration  is  the  inefficiency  that  may  result  from 
applying  constant  decisions  to  a  long  stage.  However,  generally  the  state  of  a 
system  changes  gradually,  and  thus  the  optimum  decisions  also  change  slowly. 
The  shorter  the  stage  used,  the  smaller  the  mismatch  between  the  fixed  decisions 
that  are  applied  and  the  optimum  decisions. 
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Practical  Limits  on  Short  Stages 

There  are  a  variety  of  practical  limits  on  using  short  stages.  One  universal 
difficulty  is  that  shorter  stages  impose  a  greater  computational  burden  in  solving 
the  problem.  There  may  also  be  physical  limits  of  the  system  that  prevent  rapid 
implementation  of  changes  (such  as  adjustment  of  releases  controlled  by  weirs) 
or  that  prevent  frequent  updates  of  state  variables  (such  as  efforts  to  obtain  field 
measurements  of  snow-pack  accumulation).  Though  advances  in  technology 
and  system  infrastructure  could  address  many  limitations  (such  as  using  remote 
sensing  for  snow-pack  measurements),  until  such  technology  is  in  use,  it  is 
unproductive  to  use  a  stage  that  is  shorter  than  these  limitations  allow. 

The  appropriate  stage  length  will  also  depend  on  the  characteristics  of 
system  being  modeled.  For  a  water  supply  system,  random  fluctuations  in 
demand  over  a  period  of  months  may  be  significant  but  over  a  period  of  one 
week  may  be  insignificant.  Short  term  fluctuations  become  even  less  significant 
if  there  are  significant  storage  components  which  are  not  identified  by  a  system 
model.  For  example,  a  water  distribution  system  has  storage  capacity 
downstream  of  reservoirs;  though  this  storage  may  not  be  included  in  a  system 
model,  it  still  can  contribute  to  optimal  control  by  reducing  short-term  demand 
fluctuations. 

A  special  difficulty  can  exist  when  solving  multi-objective  problems  with 
mixed  time  scales.  For  example,  the  stage  length  required  for  considering  flood 
control  objectives  may  be  inconveniently  small  for  multi-objective  problems  that 
also  consider  water  supply  or  hydropower  generation.  In  addition  to  the 
difficulty  of  finding  an  appropriate  stage  length,  it  may  also  be  difficult  to  define 
a  model  of  stream  flow  that  preserves  the  characteristics  needed  by  different 
objectives. 
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IV.  Solution  of  Lumped  Parameter  Models  Containing 
Autocorrelation  by  Gradient  Dynamic  Programming 

A.  INTRODUCTION 

Real-time  control  problems  are  a  type  of  dynamic  problem  arising  in 
many  operational  environments  such  as  process  control  and  resource 
allocation.  Often  these  problems  are  characterized  by  uncertain  inputs.  For 
such  problems,  discrete  dynamic  programming  (DDP)  is  the  only  general 
method  that  can  determine  optimal  control  solutions  without  major  and  often 
inappropriate  simplifying  assumptions.  Unfortunately,  DDP  is  limited  in  its 
ability  to  solve  complex  problems  due  to  the  exponential  growth  of 
computational  effort,  referred  to  as  the  "curse  of  dimensionality".  As  a  result, 
the  discrete  dynamic  programming  (DDP)  approach  has  been  limited  to 
solving  lumped  parameter  problems  consisting  of  only  two  or  three  state 
variables. 

Foufoula-Georgiou  and  Kitanidis  [1988]  developed  the  new  method  of 
Gradient  Dynamic  Programming  (GDP)  to  reduce  this  limitation  and  solve 
problems  consisting  of  perhaps  six  or  more  state  variables.  Though  we  can 
solve  more  complex  problems  using  GDP,  we  cannot  include  autocorrelation 
of  inputs,  such  as  stream  flow,  in  the  existing  algorithm.  I  present  an 
extension  of  the  GDP  method  that  allows  its  application  to  control  problems 
that  include  autocorrelated  inputs  or  other  information  that  can  be  used  to 
improve  estimates  of  uncertain  inputs. 

This  section  presents  justification  for  the  use  of  dynamic  programming 
and  for  the  development  of  GDP  as  an  improved  method  for  solving  dynamic 
programming  problems.  Section  2  outlines  the  basic  GDP  approach  and  Section 
3  presents  the  extended  GDP  algorithm.  The  final  section  discusses  general 
considerations  in  selecting  GDP  as  an  optimization  method. 

1.  Discrete  Dynamic  Programming 

Dynamic  programming  describes  a  class  of  optimization  techniques 
applied  to  the  control  and  planning  of  dynamic  systems  [Yakowitz,  1982]. 
Dynamic  systems  are  those  that  progressively  change  through  stages, 
representing  increments  in  either  time  or  space.  Though  frequently  confused 
with  the  more  general  concept  of  dynamic  programming,  DDP  is  but  one  class  of 
methods  available  for  solving  dynamic  programming  problems.  DDP  has  long 
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been  a  popular  optimization  technique,  especially  in  water-systems  management 
[Yakowitz,  1982;  Yeh,  1985].  The  popularity  and  success  of  this  technique  is 
primarily  due  to  its  ability  to  incorporate  nonlinear  constraints  and  objectives, 
and  its  ability  to  also  incorporate  stochastic  variables  [Willis  and  Yeh,  1987]. 

Indeed,  DDP  remains  the  only  approach  of  any  universality  in  solving  dynamic 
problems  containing  stochastic  inputs  [Yakowitz,  1982].  This  is  in  contrast  to 
deterministic  dynamic  programming  problems,  which  can  be  solved  by  a  variety 
of  approaches. 

Because  dynamic  problems  containing  stochastic  inputs  present  special 
challenges,  they  are  often  identified  as  stochastic  dynamic  programming  (SDP) 
problems.  Most  hydrologic  systems,  for  example,  present  SDP  problems  because 
they  contain  stochastic  inputs  such  as  stream  flow  and  rainfall.  Since  DDP  is  the 
only  approach  generally  applicable  to  the  solution  of  SDP  problems,  SDP  is 
commonly  synonymous  with  DDP.  However,  unlike  DDP,  SDP  does  not  specify 
a  particular  optimization  approach.  A  variety  of  approaches  have  been  used  to 
solve  specific  SDP  problems  [Karamouz  and  Vasiliadis,  1992],  though  many, 
including  GDP,  are  variants  of  DDP. 

DDP  solves  a  control  problem  by  determining  decisions  that  minimize 
expected  cost  or  other  objective  given  knowledge  of  the  current  "state"  of  a 
system.  The  "state"  of  a  system  is  a  general  term  that  describes  any  changing 
characteristics  that  influence  optimal  decisions.  Both  the  state  of  the  system  and 
optimal  decisions  evolve  from  one  stage  to  the  next  as  a  result  of  prior  decisions 
and  uncertain  inputs. 

Unlike  deterministic  systems,  a  single  control  schedule  does  not  define 
future  states  of  a  stochastic  system.  Because  future  control  decisions  depend  on 
the  uncertain  evolution  of  future  states  of  the  system,  a  stochastic  solution 
requires  the  generation  of  control  decisions  and  expected  costs  for  every  possible 
state  of  the  system  and  each  stage  [Kelman  et  ah,  1990]. 

Because  most  systems  are  characterized  by  continuous  state  variables, 
there  are  an  infinite  number  of  unique  states;  and  the  solution  consists  of  an 
infinite  number  of  control  decisions.  One  way  that  we  can  describe  an  infinite 
number  of  control  decisions  is  to  assume  that  the  expected  cost  and  control 
decisions  have  a  particular  form  of  functional  dependence  on  the  state  variables. 
For  example,  Karamouz  and  Houck  [1982]  specified  decisions  as  a  linear  function  of 
state  variables,  and  Buras  [1972]  recommended  using  a  fitted  polynomial 
function  of  high  order.  This  approach  generally  results  in  sub-optimal  solutions. 
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though  it  is  frequently  used  to  provide  a  local  approximation  in  deterministic 
solution  methods,  such  as  Differential  Dynamic  Programming  [Yakowitz,  1982; 
Jones  et  al,  1987;  Trezos  and  Yeh,  1987,  Culver  and  Shoemaker,  1993]. 

In  contrast,  DDP  methods  summarize  an  operating  policy  by  identifying 
the  expected  costs  and  optimal  control  decisions  at  each  of  a  finite  number  of 
discrete  states  and  interpolating  values  at  intermediate  states.  To  optimally 
control  a  dynamic  system  using  DDP,  an  objective  function  value  is  calculated 
for  each  feasible  transition  between  states  of  each  stage  of  a  problem.  The 
optimum  control  decisions  are  then  chosen  as  those  that  produce  the  minimum 
expected  cost. 

A  caution  in  using  DDP  is  that  one  should  be  wary  of  proposed  solutions 
until  the  deleterious  effects  of  discretization  have  been  somehow  bounded  and 
found  acceptable  [Yakowitz,  1982].  As  Gal  [1979]  pointed  out,  a  DDP 
approximation  is  optimal  only  for  the  discretized  version  of  the  model.  In  spite 
of  the  simplifying  assumptions  required  to  use  a  continuous,  global  function 
such  as  proposed  by  Buras,  a  DDP  approximation  can  perform  worse  if  the 
appropriateness  of  the  discretization  is  not  verified. 

Because  dynamic  programming  allows  us  to  avoid  specifying  in  advance 
the  functional  relationship  between  the  solution  (optimum  decisions)  and  the 
independent  variables  (state  of  the  system),  almost  any  form  of  non-linearity  in 
this  relationship  is  possible.  This  ability  is  especially  useful  in  studying  problems 
where  prior  work  has  not  provided  a  basis  for  assuming  simple  pre-determined 
functional  relationships. 

2.  Gradient  Dynamic  Programming 

Unfortunately,  DDP  problems  become  computationally  intractable  if 
many  state  variables  are  required  to  describe  the  state  of  a  system.  DDP 
describes  each  continuous  state  variable  by  a  number  of  discrete  values,  and  each 
combination  of  these  values  is  assessed  in  finding  a  solution.  As  a  result,  the 
effort  required  to  solve  a  problem  is  Tin)  Vn,  where  n  is  the  number  of 
continuous  state  variables  (referred  to  as  the  "dimension"  of  the  problem),  V  is 
some  average  number  of  discrete  values  used  to  span  each  variable,  and  Tin)  is 
the  effort  to  determine  a  solution  for  each  of  the  Vn  discrete  states.  With  each 
additional  variable  added  to  a  DDP  problem,  the  number  of  unique 
combinations  is  multiplied  by  the  number  of  discrete  values  used  to  span  the 
new  variable.  Esmaeil-Beik  and  Yu  [1984]  provide  a  good  illustration  of  this 
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difficulty  in  their  solution  of  a  reservoir  management  problem.  This  "curse  of 
dimensionality"  inherent  in  the  solution  of  DDP  problems  has  generally  limited 
its  practical  application  to  problems  described  by  only  two  or  three  state 
variables  [Yakowitz,  1982;  Johnson  et  al,  1993]. 

There  are  two  basic  approaches  available  to  reduce  the  computational 
effort  required  to  solve  a  DDP  problem.  The  first  is  to  reduce  the  number  of 
discrete  values  used  to  span  the  range  of  the  state  variables.  This  approach  may 
lead  to  a  deterioration  in  solution  accuracy,  however.  The  second  is  to  reduce  the 
number  of  variables  needed  to  describe  a  system  [Gal,  1979;  Saad  and  Turgeon, 
1988;  Saad  et  al,  1992].  For  example,  Kelman  et  al  [1990]  modeled  the  storage  of  a 
multi-reservoir  system  with  a  single  state  variable  for  storage.  However, 
aggregation  of  characteristics  can  be  inappropriate  if  we  are  to  represent  the 
system  realistically.  Beyond  these  two  general  approaches,  Johnson  et  al  [1988, 
1993]  discuss  a  variety  of  other  techniques  that  may  significantly  reduce  the 
computational  effort  of  some  problems.  These  include  solving  a  sequence  of 
problems  with  increasing  difficulty,  eliminating  from  consideration  uninteresting 
areas  of  the  state-space,  partitioning  the  original  problem  into  smaller  separable 
problems,  or  creatively  choosing  variables  used  to  model  the  system. 

To  allow  application  of  dynamic  programming  techniques  to  problems  of 
larger  dimensionality,  Kitanidis  and  Foufoula-Georgiou  [1987]  and  Foufoula- 
Georgiou  and  Kitanidis  [1988]  took  the  first  basic  approach  in  developing  GDP. 
GDP  is  able  to  significantly  reduce  the  number  of  discrete  values  required  to 
span  state  variables  while  still  achieving  a  reasonable  level  of  solution  accuracy. 
They  accomplished  this  by  using  a  more  accurate  interpolation  method  to 
determine  the  values  cost  and  decision  values  between  discrete  states.  Though 
the  effort  required  to  solve  a  problem  using  GDP  still  increases  at  an  exponential 
rate  Vn  with  dimension  n  of  the  problem,  the  base  of  the  exponent  V  can  be 
reduced.  Thus  one  can  greatly  reduce  the  number  of  unique  combinations  of 
discrete  values  that  must  be  considered. 

GDP  is  able  to  employ  a  more  accurate  interpolation  scheme  by  using  an 
interpolation  that  is  smoother  and  that  preserves  additional  information  at  each 
discrete  state.  In  addition  to  preserving  the  values  at  each  node,  GDP  also 
determines  and  preserves  the  gradients  of  these  values  with  respect  to  the  state 
variables.  These  gradients  can  be  calculated  analytically  (as  presented  here)  or 
approximated  by  a  finite  difference  approach.  To  interpolate  a  function  between 
discrete  states  while  preserving  both  the  values  and  gradients,  Kitanidis  and 
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Foufoula-Georgiou  [1987]  proposed  using  Hermite  interpolation,  developed  for 
this  purpose  by  Kitanidis  [1986].  They  demonstrated  that  with  decreases  in  the 
discretization  interval  Ax,  the  error  of  the  control  policy  and  the  cost  functions 
converge  as  (Axf  and  (Ax)4,  respectively,  using  Hermite  interpolation  versus 
Ax  and  (Axf  using  linear  interpolation. 

3.  Extension  to  Include  Autcorrelation 

A  main  purpose  of  this  work  is  to  extend  the  GDP  algorithm  to  allow  its 
application  to  problems  with  autocorrelation. 

To  model  autocorrelation,  one  frequently  includes  prior  realizations  of 
stochastic  inputs  as  state  variables.  The  probability  distribution  of  stochastic 
inputs  is  then  an  unchanging  or  "stationary"  function  of  the  state  of  the  system. 

Many  previously  solved  stochastic  control  problems  have  ignored 
autocorrelation  because  of  the  limited  number  of  state  variables  that  DDP 
methods  have  been  able  to  handle.  Control  problems  that  have  included 
autocorrelation  have  been  greatly  simplified.  Most  frequently  these 
simplifications  have  involved  the  reduction  of  significant  portions  of  a  system  to 
description  by  a  single  state  variable  [Gal,  1979;  Kelman  et  al.,  1990].  Various 
authors  have  also  applied  adaptive-control  methods  to  incorporate 
autocorrelation  while  avoiding  the  use  of  extra  state  variables  [Bras  et  al,  1983; 
Stedinger  et  al,  1984;  Karamouz  and  Vasiliadis,  1992].  Adaptive-control 
incorporates  the  effect  of  autocorrelation  by  making  transient  modifications  to 
the  probability  distributions  of  the  stochastic  inputs.  As  a  result,  adaptive- 
control  does  not  develop  global  solutions:  it  only  solves  the  optimal  control  for  a 
single  state  of  the  system  and  for  a  particular  set  of  prior  inputs.  Thus  a  new 
solution  must  be  developed  for  each  different  state,  prior  input,  or  stage. 

Compared  to  other  DDP  methods,  GDP  has  greater  flexibility  to  include 
extra  state  variables;  however,  the  existing  GDP  algorithm  is  not  general  enough 
to  describe  dependence  of  stochastic  inputs  on  state  variables  and  thus  is  unable 
to  incorporate  autocorrelation.  The  algorithm  presented  here  extends  the  GDP  to 
incorporate  this  dependence.  In  addition  to  autocorrelation,  the  extended  GDP 
algorithm  can  also  incorporate  other  correlated  information  that  conditions  the 
probability  distribution  of  stochastic  variables. 
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B.  PROBLEM  SPECIFICATION 


Dynamic  control  of  a  system  requires  that  managers  determine  control 
decisions  for  each  of  a  series  of  stages  that  define  a  planning  horizon  or  a  process. 
For  example,  real-time  management  requires  that  managers  determine  periodic 
decisions  about  how  to  operate  the  system. 

When  decisions  are  influenced  by  uncertain  inputs,  a  manager  cannot 
specify  effective  and  efficient  decisions  by  a  single  control  schedule.  Prior  values 
of  stochastic  inputs  influence  optimal  decisions,  and— without  foreknowledge— a 
manager  does  not  have  all  information  needed  to  determine  optimal  decisions 
for  a  stage  until  the  beginning  of  that  stage. 

Information  that  influences  optimal  decisions  and  that  describes  changing 
system  characteristics  specifies  the  state  of  a  system.  A  solution  of  an  optimal 
control  problem  is  the  set  of  optimal  decisions  and  expected  costs  for  each  stage 
and  state  of  a  system.  We  can  describe  these  decisions  and  costs  as  a  function  of 
the  state  of  the  system.  For  example,  optimal  water  releases  from  a  reservoir  and 
expected  risk  of  water  rationing  or  flooding  can  be  described  as  a  function  of 
current  reservoir  levels,  stream  flows,  and  season. 

1.  Modeling  a  Control  Problem 

We  can  measure  the  performance  of  decisions  with  an  objective  function 
that,  in  its  most  general  form,  can  be  represented  as  any  function 

f  =  /(x,u,w,y) 

where  x  is  a  vector  of  independent  state  variables  describing  the  system  at  the 
beginning  of  the  current  stage,  u  is  a  vector  of  decision  variables  applied  during 
the  stage,  w  is  a  vector  of  random  variables  whose  values  are  realized  during 
the  stage,  and  y  is  a  vector  of  state  variables  describing  the  system  at  the  end  of 
the  stage.  Vectors  x  and  y  describe  the  same  state  variables  in  adjacent  stages; 
i.e.,  y,  =  x,+i  in  adjacent  stages  t  and  r+1 . 

Variables,  u,  w,and  y  are  implicit  functions  of  state  variables,  x. 
Autocorrelation  is  modeled  by  making  the  distribution  of  random  variables 
dependent  upon  the  state  of  the  system,  thus  w(x) .  A  transition  function, 
y  =  T(x,u,w)  ,  models  the  evolution  of  a  system  from  state  x  to  state  y.  The 
control  solution,  u*,  is  specified  as  a  function  of  x,  thus  u*(x). 
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In  order  to  incorporate  stochastic  variables  in  a  model,  only  a  limited 
number  of  parameters  can  be  used  to  model  a  probability  distribution.  We  can 
parameterize  a  distribution  either  by  using  its  moments  or  by  using  discrete 
probability-weighted  realizations  that  span  the  distribution.  In  this  work,  I 
estimate  w  by  discrete  realizations  w(fc)  with  probability  pk.  Likewise,  because 

y  depends  on  w,  its  values  are  uncertain  and  estimated  by 

y(jk)  =  T(x,u,W(jt))  . 

using  the  same  probability  weights. 

Foufoula-Georgiou  and  Kitanidis  [1988]  developed  the  GDP  method  for 
stage-wise  decomposable  problems  with  the  more  limited  objective  function 

f  =  C(x,u)  +  F(y) 

that  expresses  the  objective  as  a  total  "cost"  composed  of  a  current  deterministic 
cost  C  and  an  expected  future  cost  F.  F  is  only  a  function  of  y  since  the  state 
of  a  system  specifies  all  information  needed  to  determine  expected  costs  of  future 
decisions.  For  stochastic  problems,  y  and  F(y)  are  uncertain  and  the  objective 
is 


/  =  C(x,u)  +  X  {PkF(y(k))}  • 

k 

2.  Solving  a  Control  Problem 

We  can  determine  an  optimum  control  solution  by  recursively  solving  the 
system 


F(t)(  y)  =  /(f+i)*(x) 


where 


/(,+D*(x)  =  nun„{  C(,+i)(x,u)  +  F(/+i)(y) }. 

is  from  the  solution  of  the  control  problem  in  the  previous  recursion.  The 
previous  recursion  solves  the  optimal  decisions  for  the  following  stage,  and  the 
solution  of  the  problem  progresses  backwards  through  time  or  space.  The  last 
stage  is  solved  by  the  first  recursion  using  an  assumed  function  F^(y) 
representing  costs  expected  beyond  the  end  of  the  planning  horizon. 
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If  the  planning  horizon  is  very  long  or  if  future  costs  are  discounted,  a 
specified  cost  function  F^(y)  may  have  a  negligible  impact  on  current  decisions 
and  could  be  zero  or  any  arbitrary  value.  On  the  other  hand,  if  the  cost  function 
%)(y)  significantly  affects  current  decisions,  it  must  have  reasonable  values.  In 
practice,  F(W)( y)  commonly  specifies  expected  costs  that  are  the  same  for  all  final 
states;  in  other  words,  the  state  of  the  system  at  the  end  of  the  planning  horizon 
doesn’t  matter.  Less  commonly,  F(^( y)  specifies  costs  that  vary  with  the  final 
state.  For  example,  if  planned  maintenance  requires  a  certain  final  state,  F(#)(y) 
could  specify  a  high  cost  for  any  final  state  other  than  that  required.  In  their 
study  of  a  hydropower  system,  Reiman  et  al.  [1990]  used  an  assumed  cost 
function  based  on  the  potential  energy  of  water  remaining  in  reservoirs. 

In  some  well  studied  problems,  we  may  know  the  general  form  of  the 
function  F(y).  For  example,  if  we  can  assume  that  expected  costs  are  a  linear 
function  of  state  variables,  then  F(y)  is  defined  by  a  few  parameters  that  define 
the  linear  relationships.  In  many  problems,  however,  we  do  not  know  a  general 
form  of  F(y).  In  these  problems,  we  can  use  a  DDP  method  that  approximates 
F(y)  by  discrete  values  spanning  the  state  space.  For  each  discrete  state  x(t), 
DDP  solves 


/*(x(i-))  =  min„{  C(x(/),u)+  £  {p*F(y  (*,,•))}  } 

k 

where  y^,-)  =  T(X(,),u,W(^,))  ,and  w^,)  =  W(*)(X(,))  .  Once  we  have  determined 
optimal  decisions  u*(x(l))  and  values  of  the  objective  /*(x(i))  for  each  x(i), 
interpolation  between  these  values  create  piece-wise  continuous  functions  for 
u*(x)  and  /*(x) . 
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C  THE  GDP  METHOD 


1.  Notation 

Vectors  are  column  vectors.  The  transpose  of  a  column  vector  is  a  row 
vector.  Derivatives  of  a  vector  with  respect  to  a  scalar  is  a  vector  with  the  same 
dimension.  The  derivative  of  a  scalar  function  with  respect  to  a  column  vector  is 
a  row  vector.  Thus, 


but 


Ml(x)' 

du\/dx 

U(x)  = 

U2(X) 

• 

t/U 

^  3*" 

du2/dx 

9 

_  ft  _ 

Ax)  ^  —  f*  -  \fx\>  fxi>  ] 


where  u  is  a  vector  function  of  scalar  x  and  f  is  a  scalar  function  of  vector  x. 

It  follows  from  these  conventions  that  the  derivatives  of  a  m-dimensional 
vector  function  with  respect  to  a  ^-dimensional  vector  is  a  (m  x  n)  matrix  of 
derivatives  where  the  dependent  vector  function  is  indexed  along  rows  and  the 
independent  vector  indexed  along  columns.  Second  derivatives  of  scalar 
functions  with  respect  to  vectors  are  also  matrices  with  the  first  vector  indexed 
along  rows  and  the  second  vector  indexed  along  columns.  Thus 


3u 

dx 


dui/dxi  dui/dx2 
durfdxi  du2/dx2 


ouox 


are  matrices  of  derivatives,  respectively,  of  vector  function  u(x)  with  respect  to 
vector  x  and  of  scalar  function  f  with  respect  to  vectors  x  and  u. 

When  /  =  /(x,u(x))  ,  the  definition  of  the  total  derivative  with  respect  to 

x  is 


V/= 


dx  dx  du  3x 


I  distinguish  total  and  partial  derivatives  of  scalar  functions,  respectively,  by  the 
use  of  the  subscripted  symbol  V  versus  subscripts  directly  on  the  function.  In 
contrast,  total  and  partial  derivatives  of  vector  functions  are  distinguished  by  the 
use  of 'd'  and 'd'  respectively,  and  will  not  use  the  symbol  V  or  subscripts. 
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2.  Method 

The  following  development  of  the  GDP  method  parallels  that  of  Foufoula- 
Georgiou  and  Kitanidis  [1988]  with  the  exception  that  I  use  index  k  to  indicate  a 
discrete  realization  of  stochastic  variables  and  index  t  to  indicate  a  stage. 

The  goal  is  to  provide  rules  that  allow  the  system  operator  to  make  real¬ 
time  operating  decisions.  These  rules  determine  optimal  decisions  u*  that 
minimize  an  objective  function  f ;  in  other  words 

min„{/(x,u,w,y) } 

subject  to  the  set  of  active  constraints 

Gu  =  d(x) 

where  G  is  a  (p  x  m)  scalar  matrix  of  rank  p,  and  d(x)  is  a  vector  of  the  m  right 
hand  sides.  The  general  form  of  an  inequality  constraint  is 

gTu  <  d(x) 

where  g  is  a  vector  of  coefficients  and  d  can  be  a  function  of  x.  When  binding, 
an  inequality  constraint  becomes  an  equality  constraint  in  the  active  set. 

Using  the  method  of  Lagrange  multipliers,  optimal  decisions  that  respect 
the  constraints  are  also  those  that  minimize 

h(u,X)  =f+  A.T(Gu  -  d) 

with  respect  to  u  and  X,  where  X  is  a  vector  of  Lagrange  multipliers.  Taking 
derivatives  of  the  above  equation,  optimal  controls  are  defined  by  the 
simultaneous  equations 

(Vu/i)T  =  (Vu/)T  +  XJG  =  0 
(Vx/i)T  =  Gu  -  d  =  0  . 


Though  these  equations  can  verify  optimal  decisions  u*,  they  cannot  be 
used  to  directly  solve  u*.  Objective  function  f  is  derived  from  the  piece- wise 
continuous  function  F,  and  we  do  not  know  which  "piece"  of  this  function  to  use 
until  after  specifying  the  final  state  y  resulting  in  part  from  control  decisions. 
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Since  evaluation  of  the  above  equations  requires  this  knowledge,  we  cannot,  in 
turn,  determine  optimal  decisions  until  the  appropriate  "piece"  of  the  expected 
cost  function  is  specified. 

Instead,  given  an  initial  set  of  control  decisions,  the  algorithm  determines 
incremental  changes  Au  that  iteratively  approach  u*  by  applying  a  Newton- 
type  method.  Subsequently,  by  determining  u*  for  each  state  X(/),wecan 
develop  the  implicit  control  function  u*(x)  and  cost  function  /*(x) .  Using 
/*(x)  as  the  prior  stage's  expected  cost  function,  F(y),  then  allows  us  to 
determine  a  solution  for  the  prior  stage.  This  new  solution  then  represents  a 
planning  horizon  one  stage  longer. 

3.  Algorithm 

The  algorithm  consists  of  three  nested  cycles,  the  most  basic  of  which  is 
the  determination  of  u*.  The  optimal  decisions  u*  apply  only  to  a  discrete 
initial  state  x©  during  recursion  t.  For  each  additional  recursion,  the  algorithm 
solves  u*  for  every  x(()  until  it  has  performed  enough  recursions  to  span  a 
planning  horizon  of  N  stages  or  to  allow  control  decisions  to  converge.  Figure 
IV-1  illustrates  the  three  cycles  and  the  steps  of  the  algorithm. 

Step  1:  Discretize  the  State  Space 

A  control  solution  specifies  u*(x)  and  /*(x)  as  functions  of  state 
variables  x  =  (jti, ...,  xn)T  .  Since  we  generally  do  not  know  the  form  of  these 
functions,  the  algorithm  solves  piece-wise  continuous  functions  by  calculating 
values  at  a  number  of  discrete  states  x(j)  and  approximating  intermediate  values 
by  interpolation. 

The  accuracy  and  efficiency  of  the  GDP  solution  depends  on  the  selection 
of  an  appropriate  number  of  discrete  values  X(,) .  Selection  of  discrete  states 
involves  trade-offs  between  solution  accuracy  and  computational  effort  [Esmaeil- 
Beik  and  Yu,  1984].  Efficient  placement  of  the  x©  also  requires  judgment:  state 
variables  that  induce  rapid  changes  in  a  function  need  finer  discretization;  state 
variables  that  have  a  consistent  effect  on  a  function  allow  coarser  discretization. 
We  should  verify  a  discretization's  appropriateness  prior  to  accepting  a  solution 
[Yakowitz,  1982]  such  as  by  observing  the  effect  of  using  different  discretizations. 

Step  2:  Specify  an  Expected  Cost  Function 

The  expected  cost  function  F(y)  is  a  piece-wise  interpolated  function  with 
known  value  F  and  gradients  VyF  at  each  discrete  state  y (,>  For  the  first 


81 


recursion,  we  must  provide  values  and  gradients  that  describe  F(/v>( y),  the 
function  for  costs  expected  beyond  the  end  of  a  planning  horizon  of  N  stages. 
For  remaining  recursions  t =N-\,N-2, ...,  1  ,  the  algorithm  provides  values  and 
gradients  from  the  previously  solved  recursion  where  y (,),<  =  X(o,r+i  and 

Given  a  discretization  and  values  at  each  discrete  state,  GDP  is  able  to 
interpolate  values  at  intermediate  points  using  Hermite  interpolation(Appendix 
B).  Hermite  interpolation  preserves  both  values  and  gradients  at  the  discrete 
states  that  bound  intervals  [Kitanidis,  1986].  Foufoula-Georgiou  and  Kitanidis  [1988, 
Appendix  A]  outline  the  Hermite  interpolation  method  in  greater  detail,  with 
correction  that  their  equation  (B 7)  should  be  T  =  (1  -  r\j)(  1  -  3 T)})  for  s  =;. 

Step  3:  Select  a  Discrete  State  and  Initial  Control  Decisions 

In  order  to  develop  functions  u*(x)  and  /*(x) ,  the  algorithm  determines 
optimal  values  for  u*,  f* ,  and  for  each  discrete  state  x©  by  steps  4  to  9. 
Given  an  average  number  V  of  discrete  values  spanning  each  of  the  n  state 
variables,  the  total  number  of  discrete  states  is  approximately  Vn.  We  may  be 
able  to  reduce  this  number  and  achieve  greater  solution  efficiency  by  eliminating 
from  consideration  infeasible  combinations  of  state  variable  values  [Johnson  et  al, 
1988, 1993]. 

With  the  specification  of  each  discrete  state  X(,),  we  can  select  an  initial 
control  decision.  By  selecting  initial  decisions  that  are  close  to  the  optimum,  we 
can  improve  the  search  efficiency  and  likelihood  of  convergence.  Also,  if  the 
nature  of  the  problem  suggests  initial  decisions  that  satisfy  all  constraints,  these 
decisions  will  be  feasible  and  step  4  can  be  skipped.  If  no  feasible  decisions  are 
apparent,  step  4  will  test  a  set  of  arbitrary  decisions  for  feasibility  prior  to 
independently  determining  initial  feasible  decisions. 

Step  4:  Determine  Initial  Feasible  Decisions 

If  the  initial  decisions  specified  in  step  3  violate  any  of  the  constraints,  the 
"phase  1"  procedure  of  linear  programming  can  determine  feasible  control 
decisions.  Luenberger  [1984]  describes  the  phase-1  procedure. 

Step  5:  Calculate  Improvement  in  the  Decisions 

The  iterative  improvement  Au  is  the  solution  to  the  system  of  equations 
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where 


V,(v,/)T 

gt 

Au 

'-v/ 

G 

0 

.  X  . 

0 

Vuf  =  c„  +  X  ( pf  ?.<« ) 

V»(V/|T  =  cu„  +  t  mP„m )  ^  • 

Matrix  G  is  from  the  initial  active  set  of  step  4  or  the  updated  set  of  step  6  or  7. 
Vector  A.  contains  the  p  Lagrange  multipliers  needed  to  check  the  active 
constraint  set  in  step  7.  Foufoula-Georgiou  and  Kitanidis  [1988]  adapted  this 
approach  from  the  active  set  method  outlined  by  Luenberger  [1984,  chapter  11]. 
The  new  solution  is  u'  =  uo  +  Au  where  uo  is  the  vector  of  old  decisions. 

Step  6:  Check  Feasibility  of  Decisions 

If  an  improvement  Au  makes  new  decisions  infeasible,  the  active 
constraint  set  is  augmented  with  the  violated  constraint.  The  algorithm  then 
identifies  new  feasible  decisions  u’  by  projecting  uo  on  the  domain  defined  by 
the  updated  active  constraint  set.  The  new  control  decisions  are  u’  =  uo  +  5u 
where 

8u  =  GT(GGT)-1(d  -  Gu0) 

[Foufoula-Georgiou  and  Kitanidis,  1988,  Appendix  C]. 

An  empty  feasible  decision  space  may  result  if  the  algorithm  adds  more 
than  one  constraint  to  the  active  set  at  a  time.  When  an  improvement  Au 
violates  more  than  one  constraint,  the  algorithm  adds  to  the  active  set  only  the 
first  constraint  encountered  when  moving  from  uo  in  the  direction  Au.  The 
distance  from  uo  to  each  new  constraint  is 

y  =  (d-  gTu)  /  gTAu 

and  the  first  violated  constraint  is  that  with  the  smallest  distance  y . 

Step  7:  Check  if  Optimum  Decisions  are  Found 

Control  decisions  are  close  to  optimum  when  the  working  set  of 
constraints  is  still  correct  for  the  updated  decisions  and  changes  are  small 
(Jau|  <  8  or  AuTAu  <  8).  If  either  condition  is  not  met,  the  algorithm  returns  to 
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step  5  after  removing  any  unnecessary  constraints.  Unnecessary  constraints  are 
those  with  a  negative  Lagrange  multiplier  (A  <  0). 

Step  8:  Calculate  Values  of  Total  Cost  Function  and  Derivatives 

Given  the  optimal  decisions  u*,  the  algorithm  can  now  calculate  the 
optimum  value  J*(x)  =  min„{/(x,u,w,y) }  and  gradients  V*/*  for  the  current 
initial  state  x®.  These  values  are  constants  for  the  current  xg)  since  the  optimal 
decisions  are  also  constants  for  X(j) . 

The  optimal  expected  value  of  the  objective  function  is 

/*(x(0)  =  C(x(l>u*)+  X  IpiNjm*)  > 

k 


where  y&o*  =  T(x(l),u*,w(^i))  . 

Calculation  of  the  matrix  of  derivatives  requires  du*/dx.  The 

algorithm  calculates  the  matrix  of  values  du*/dx  by  solving  the  system  of 
equations: 


Vu(VuOT 

Gt 

du*/dx 

-D 

G 

0 

.  dX./dx  . 

Vxd 

i  n  _  ,  a.vTv  .  r  /Sy  ,  3w<«  \ , 

where  D  =  C»*  +  f  < PtF„m  (5^  +  5^-  -5-J ) 

and  Vxd  =  dx  • 


Appendix  C  provides  the  derivation  of  the  above  system  of  equations.  Given 
du*/dx. 


W)T  =  <v«oT^ 


+  Cv 


jyw  ^w(^)  \  1 

dw(k)  dx  I 


Step  9:  Completing  the  Solution  and  Results 

Given  the  optimal  controls  u*,  the  algorithm  returns  to  step  3  to  repeat 
steps  3  to  8  for  each  initial  state  X(/>.  Once  the  optimal  values  of  u*,  du*/dx,  f*  , 
and  Vxf1  have  been  found  for  each  discrete  x^,  the  algorithm  can  specify  u*(x) 
and  J*(x)  for  any  x  using  Hermite  interpolation. 

The  piece-wise  continuous  functions  u*(x)  and  /*(x)  are  the  solution  to 
a  control  problem  having  a  planning  horizon  that  is  now  one  stage  longer  than 
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that  of  the  previous  recursion.  Given  u*(x)  and  /*(x) ,  the  algorithm  returns  to 
step  2  to  repeat  steps  2  to  8  for  each  remaining  recursion  t  =  N-l,  N-2, ...  1  . 

Unlike  many  optimal  control  algorithms  that  require  a  forward  recursion  to  solve 
for  the  optimal  controls  from  the  expected  cost  function,  the  backwards  recursion 
generates  both  u*(x)  and  /*(x)  for  the  entire  planning  horizon;  therefore,  the 
solution  is  complete  following  the  backwards  recursion. 
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Fig.  IV-1.  Flow  chart  for  the  GDP  algorithm 
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D.  COMMENTS  AND  CONCLUSIONS 

To  solve  a  control  problem  using  simulation  or  optimization  requires  that 
we  develop  a  simplified  model  of  a  system.  GDP  makes  it  possible  to  solve  DDP 
problems  using  models  that  are  less  simplistic  than  past  efforts.  Though  this 
makes  GDP  an  appropriate  optimization  method  for  many  problems,  it  is  not 
necessarily  the  best  DDP  technique  for  all  problems.  GDP  is  complicated  to 
implement,  and  becomes  useful  only  when  it  is  not  possible  to  solve  a  problem 
by  other  techniques  without  excessive  simplification.  Given  the  limitations  of 
most  DDP  methods  however,  this  occurs  frequently. 

There  are  some  specific  characteristics  of  GDP  that  we  should  consider 
prior  to  applying  it  to  a  problem.  In  addition  to  its  added  complexity,  GDP 
should  also  be  evaluated  in  terms  of  its  efficiency  and  ability  to  converge  on  a 
solution.  I  discuss  these  briefly  prior  to  presenting  my  conclusions. 

1.  Efficiency 

Johnson  et  al.  [1988, 1993]  compared  the  efficiency  of  various  spline 
interpolations  against  the  Hermite  interpolation  used  in  GDP  and  against  linear 
interpolation.  Spline  techniques  of  Johnson  et  al. ,  like  the  Hermite  interpolation, 
seek  to  improve  accuracy  and  smoothness  of  estimates;  however,  they  differ  in 
that  they  do  not  use  additional  information  from  gradients. 

Both  spline  and  Hermite  interpolations  demonstrated  the  ability  to  solve  a 
4-dimensional  problem  with  less  than  1%  the  effort  and  equivalent  accuracy  of  a 
linear  method.  They  also  demonstrated  that  Hermite  interpolation  is  more 
accurate  than  spline  interpolation  for  the  same  discretization  of  the  state  space. 
On  the  other  hand,  Hermite  interpolation  requires  extra  computational  effort  to 
obtain  gradients;  and  they  determined  that,  for  their  problem,  GDP  was 
marginally  less  efficient  than  their  spline  interpolation  method. 

While  both  the  Hermite  and  spline  interpolations  can  significantly 
improve  the  accuracy  and  efficiency  of  a  DDP  approach,  it  is  not  clear  which  of 
these  is  better  for  a  particular  application.  Though  it  is  not  a  purpose  of  this 
work  to  compare  the  performance  of  these  different  approaches,  a  few 
observations  are  in  order. 

In  general,  the  smoother  the  function  and  the  greater  number  of  discrete 
values  used  to  span  the  state-space,  the  less  valuable  becomes  information 
provided  by  gradients.  I  expect  that  Hermite  interpolation  is  especially 
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appropriate  for  irregular  functions  (though  stochastic  problems  are  generally  less 
irregular  due  to  averaging  process  involved  in  determining  expectations).  I  also 
expect  that  Hermite  interpolation  is  more  appropriate  for  functions  that  must  be 
spanned  by  a  few  discrete  states,  as  may  be  required  when  the  number  of  state 
variables  is  large  or  when  computational  resources  are  limited.  In  an  extreme 
case  using  just  two  discrete  values  to  span  the  range  of  a  state  variable,  Hermite 
interpolation  should  generally  provide  a  better  approximation  of  the  function. 

Though  the  accuracy  of  an  interpolated  function  is  improved  more  by  the 
addition  of  a  discrete  state  than  by  knowledge  of  a  gradient,  more  gradients  are 
generally  obtained  for  an  equivalent  amount  of  computational  effort.  Whether 
calculating  gradients  or  calculating  adding  additional  discrete  states  is  more 
efficient  may  be  difficult  to  know  in  advance. 

2.  Convergence 

Past  efforts  to  approximate  an  expected  cost  function  by  analytic  functions 
have  often  relied  on  polynomials  of  high  order.  However,  fitting  high-degree 
polynomials  to  data  over  the  entire  state  space  can  induce  severe  oscillations  in 
the  fitted  function.  Optimization  applied  to  such  a  fitted  function  may 
incorrectly  seek  out  local  optima  induced  by  these  oscillations  [Johnson  et  al, 
1993].  Unlike  a  polynomial  fit,  DDP  avoids  this  problem  since  it  does  not  try  to 
fit  any  predetermined  functional  form  to  data.  Instead  it  accepts  values  at  the 
discrete  states  as  precise  and  applies  interpolation  to  specify  intermediate  values. 

However,  DDP  still  must  rely  on  an  analytic  function  to  interpolate 
between  discrete  states.  The  most  common  functions  are  the  zero'th-order 
nearest  neighbor  and  the  first-order  linear  interpolations.  Hermite  and  spline 
techniques  are  simply  more  sophisticated  approaches.  GDP  is  a  DDP  method 
that  uses  the  Hermite  interpolation  to  determine  third  or  higher-order 
polynomials  that  preserve  the  function  values  and  gradients  at  endpoints  of  a 
hyper-cube  defined  by  discrete  states.  Unfortunately,  even  third-order 
polynomials  can  create  local  minima  that  may  incorrectly  be  sought  out  by  an 
optimization  routine.  Because  of  the  order  of  the  polynomial,  at  least  one 
artificial  oscillation  is  possible.  Though  oscillations  will  generally  have  low 
amplitudes,  recursive  application  of  an  optimization  routine  can  amplify  them. 

If  an  optimization  effort  fails  to  converge  due  to  oscillations  of  the 
Hermite  interpolation,  there  are  a  variety  of  approaches  that  we  can  use  to 
achieve  convergence.  The  simplest  is,  if  possible,  to  choose  an  initial  solution 
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that  is  closer  to  optimality.  This  may  prevent  or  limit  the  growth  of  artificial  local 
optima  induced  by  oscillations.  Another  approach  is  to  use  known 
characteristics  of  the  solution  to  constrain  either  discrete  or  interpolated  values. 
For  example,  if  an  expected  cost  function  changes  monotonically  over  certain 
intervals,  a  problem  can  be  constrained  to  specify  that  no  values  interpolated  in  a 
hyper-cube  lie  outside  the  range  defined  by  bounds  of  that  hyper-cube. 

3.  Conclusions 

Incorporating  autocorrelation  into  GDP  allows  us  to  solve  more  realistic 
and  practical  problems  than  previously  possible  by  other  DDP  methods.  Because 
GDP  permits  us  to  develop  global  solutions  that  specify  optimal  control 
decisions  for  every  possible  state  of  a  system,  we  have  a  greater  opportunity  to 
observe  general  solution  characteristics.  This  may  allow  us  to  develop  general 
guidelines  for  control  of  other  problems  and  to  justify  convenient  simplifying 
assumptions  in  developing  models. 

The  extended  GDP  technique  presented  here  can  be  used  to  incorporate 
any  information  that  conditions  the  probability  distribution  of  stochastic  inputs. 
In  addition  to  autocorrelation,  there  is  often  other  correlated  information  that  we 
can  use  to  improve  predictions  of  future  inputs.  This  information  is  important 
for  determining  optimal  control  decisions  and  should  be  included  in  models 
used  to  determine  optimal  control.  The  extended  GDP  algorithm  presented  here 
is  general  and  allows  addition  of  any  state  variable  to  a  model  that  uses  them  to 
condition  the  distribution  of  stochastic  variables. 
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V.  Discussion  and  Concluding  Remarks 

A.  EXTENSIONS  OF  THIS  ANALYSIS 

The  process  of  replicating  a  system  by  an  abstract  model  involves 
significant  simplification.  In  the  model  presented  here,  I  have  used  several 
simplifications  that  may  be  inappropriate  for  the  EBMUD  system. 

Improved  Modeling  of  Physical  Characteristics 

Snow  accumulates  in  the  Sierra  Nevada  Mountains  during  wet  winter 
months  and  subsequently  feeds  stream  flow  during  the  spring  and  early 
summer.  As  a  result,  measurements  of  snowpack  accumulation  can  greatly 
improve  estimates  of  stream  flow  and  should  be  included  in  a  model  of  the 
EBMUD  system  used  to  determine  optimal  control.  Also,  the  system  consists  of 
two  reservoirs  in  the  Sierras  and  several  small  terminal  reservoirs.  While  it  may 
be  reasonable  to  lump  the  storage  of  some  of  these  reservoirs  in  the  same  state 
variable,  the  appropriateness  of  these  simplifying  assumptions  should  be 
checked.  Similarly,  the  groundwater  component  of  the  proposed  system  has 
been  reduced  to  a  "bathtub"  model  that  may  not  be  appropriate  for  this 
conjunctive-use  system. 

Addition  of  Economic  Costs 

In  addition,  I  have  not  included  estimates  of  actual  economic  costs  arising 
from  water  rationing.  Economic  costs  become  important  when  balancing  the 
costs  and  benefits  of  different  objectives.  Even  with  only  a  water  supply 
objective,  costs  of  rationing  must  be  balanced  against  costs  of  groundwater 
pumping.  To  provide  an  appropriate  basis  for  including  different  criteria, 
reasonable  measures  of  economic  costs  can  be  developed. 

Effects  of  Discretization 

There  are  various  aspects  of  dynamic  programming  that  we  should 
consider  before  accepting  a  solution.  As  pointed  out  by  Yakowitz  [1982],  "One 
should  be  wary  of  proposed  solutions  until  the  deleterious  effects  of 
discretization  have  been  somehow  bounded  and  found  acceptable."  While  GDP 
greatly  increases  the  solution  accuracy  for  a  particular  discretization, 
discretization  still  adds  some  amount  of  error  to  the  solution.  I  should  evaluate 
this  error  to  ensure  an  appropriate  discretization  for  the  solution. 
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As  pointed  out  by  Johnson  et  al.  [1993],  if  the  general  form  of  the  expected 
cost  function  or  other  characteristics  of  the  model  solution  is  known,  this 
knowledge  may  be  used  to  select  appropriate  discrete  values  of  state  variables. 
Thus,  it  may  also  be  possible  to  provide  some  general  guidance  on  the  selection 
of  appropriate  discrete  values. 

Other  Modeling  Approaches 

Finally,  the  conjunctive  use  problem  may  be  solvable  by  simpler  methods 
than  presented  here,  thus  allowing  the  development  of  more  realistic  models.  As 
discussed  earlier,  we  may  be  able  to  identify  appropriate  functions  that  can  be 
used  to  avoid  discretization  of  some  or  all  state  variables.  Alternatively,  we  may 
be  able  to  break  the  problem  into  simpler  sub-problems  or  restrict  our  solution  to 
avoid  considering  uninteresting  states  or  decisions. 

B.  CONCLUSIONS 

This  work  has  expanded  the  ability  of  dynamic  programming  techniques 
to  solve  a  wider  variety  of  stochastic  problems  by  Gradient  Dynamic 
Programming.  By  the  expanded  algorithm  presented  here,  I  have  addressed  a 
problem  that  contains  autocorrelation  without  resorting  to  heuristic  constraints 
on  the  control  of  groundwater. 

The  solution  of  a  generic  conjunctive  use  problem  has  allowed  us  to 
develop  some  initial  insight  into  the  optimal  control  of  these  problems.  This 
problem  will  form  the  basis  for  problems  that  model  the  system  more 
realistically.  This  will  promote  the  goal  of  providing  practical  guidance  to 
planners  and  managers  of  conjunctive  use  systems  attempting  to  use  a 
conjunctive-use  approach. 

The  Gradient  Dynamic  Programming  method  presented  here  also  has 
application  beyond  that  of  conjunctive  use.  It  can  be  a  useful  tool  for  any 
stochastic  dynamic  programming  problem,  even  beyond  those  of  water  resource 
management. 
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C  DEFINITION  OF  VARIABLES 


u  vector  of  control  decision  variables  up  j=l,...n 

u*  vector  of  optimal  control  decisions 

x  vector  of  state  variables  Xj,  j  =  1  ,...n 

X(,-)  values  of  state  variables  at  discrete  state  (i),  i  =  1....V  " 

y  vector  of  state  variables  at  the  end  of  the  current  decision  period 

w  vector  of  stochastic  variables  Wj,  j  =  1  ,...n 

W(jt)  values  of  stochastic  variables  at  discrete  realization  ( k ) 

T  vector  of  transition  functions  such  that  y=T(x,u,w) 

g,  g  function  and  coefficient  vector  defining  the  left-hand  side  of  constraints 
G  array  of  left-hand  side  coefficients  of  currently  active  constraint  set 
d  function  defining  the  right-hand  side  of  constraints 
d  vector  of  right-hand  side  coefficients  of  currently  active  constraint  set 
f  total  cost  function,  /(x,u,w)  =  C  +  F 
f*  optimal  total  cost  function,  j*  =/(x,u*,w) 

C  current  cost  function,  C(x,u) 

F  expected  future  cost  or  "cost-to-go"  function,  F(y) 

Ew{ }  expected  value  operator  with  respect  to  w 
Pk  probability  weight  for  W(*) 

t  index  for  stage 

t  index  for  month 

wt  a  single  stochastic  variable  w  at  time  t 

(Ox  log-transformed  stream  flow  for  month  t  using  the  3-parameter  model 
mean  of  (Ox 

Gx  standard  deviation  of  Oh. 

Xx  third  parameter  in  the  three-parameter  lognormal  distribution 
Px  lag-one  autocorrelation  coefficient  of  log-transformed  stream  flows 
£  random  variable  with  normal  (0,1)  distribution 


q  quantile  of  a  probability  distribution 

d  partial  derivative  operator 

d  total  derivative  operator  (defined  by  the  chain-rule) 
fx  shorthand  notation  for  partial  derivatives  (here  of  f  with  respect  to  x) 
used  with  scalar  functions  (eg.  f ,  C,  F) 


92 


V/  shorthand  notation  for  total  derivatives  (here  of  f  with  respect  to  x) 
used  with  scalar  functions 
%  partial  autocorrelation  coefficient 

c  standardization  constant  applied  to  expected  total  cost  function 

N  number  of  stages  in  the  planning  horizon 

R  matrix  of  partial  derivatives  dw/9x 

rjf  elements  of  R 

P  matrix  of  correlation  coefficients 

a  discount  rate 

P  discount  factor 

V  average  number  of  discrete  values  to  span  each  state  variable 
n  dimension  of  problem  (i.e.,  number  of  state  variables) 

T  effort  to  determine  u*  for  one  discrete  state  X(q  and  stage  t 
Ax  discretization  interval 

h  total  cost  function  augmented  with  constraints  by  Lagrange  multiplier 
method 

A  Lagrange  multipliers  that  minimize  h 

D  shorthand  notation  used  for  a  defined  matrix  calculation 
y  distance  from  solution  uo  to  the  first  new  constraint  in  direction  Au 
4>  Hermite  weights  applied  to  function  values 

y/  Hermite  weights  applied  to  function  derivatives 

xa,  xb  upper  and  lower  hypercube  bounds 

distance  between  x  and  xb  in  dimension  k 
i)k  distance  between  x  and  X(t)  in  dimension  k 
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Appendix  A 

Constraining  the  Optimization  Problem 

The  goal  is  to  provide  rules  that  allow  the  system  operator  to  make  real¬ 
time  operating  decisions.  These  rules  determine  the  control  decisions  u*  that 
minimize  an  objective  function  f ;  in  other  words 

minu(/(x,u,w,y) }  subject  to  active  constraints  Gu  =  d(x) 

where  G  is  a  (p  xm)  matrix  of  rank  p,  and  d(x)  is  a  vector  of  the  m  right  hand 
sides.  The  general  form  of  each  inequality  constraint  is 

gTu  <  d(x) 

where  g  is  a  vector  of  coefficients  and  d  can  be  a  function  of  x.  When  binding, 
an  inequality  constraint  becomes  an  equality  constraint  in  the  active  set. 

An  empty  feasible  decision  space  may  result  if  the  algorithm  adds  more 
than  one  constraint  to  the  active  set  at  a  time.  When  an  improvement  Au 
violates  more  than  one  constraint,  the  algorithm  adds  to  the  active  set  only  the 
first  constraint  encountered  when  moving  from  uo  in  the  direction  Au.  The 
distance  from  uo  to  each  new  constraint  is 

7  =  (d  -  gTu)  /  gTAu 

and  the  first  violated  constraint  is  that  with  the  smallest  distance  7. 

1.  Formulation  of  the  Constraints 

A  constraint  on  the  decision  variables  u  has  the  general  form 

g(u)  <  d(x,xv) 

where  g(u)  is  a  function  of  u  and  right  hand  side  d  is  a  function  of  x  and  w. 
Function  d(x, w)  is  not  dependent  on  y  since  y  can  be  replace  by  the  transition 
function  y  =  T(x,u,w  ).  The  "less-than-or-equal-to"  form  is  standard  for  non¬ 
linear  optimization  and  allows  consistent  interpretation  of  the  Lagrange 
multipliers  when  verifying  the  active  constraint  set.  Equality  constraints  are 
always  active. 
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Assuming  a  system  with  linear  dynamics  (i.e.,  no  second  order  terms  in 
the  transfer  function),  the  terms  containing  u  can  be  separated  from  those  with 
x  and  w.  Thus,  the  constraints  can  be  represented  by 

gTu  <  d 

where  g  is  a  vector  of  coefficients.  Limited  non-linear  dynamics  can  be 
represented  by  this  representation  of  constraints  as  long  as  non-linear  terms  are 
limited  to  d(\, w) . 

The  active  set  method  used  to  solve  the  problem  relies  on  continuously 
updating  a  set  of  constraints  limited  to  those  that  bind  the  current  solution. 
Constraints,  when  active,  become  equality  constraints,  and  the  entire  active  set 
can  be  represented  by 


Gu  =  d 

where  G  is  a  scalar  matrix  of  coefficients,  and  d  is  a  vector  of  right  hand. 

Vector  d  is,  most  generally,  a  function  of  x  and  w .  The  dependency  on 
w  can  be  removed  by  assuming  an  appropriately  conservative  value  of  w.  By 
over-constraining  u,  a  solution  will  always  (or  practically  always)  be  feasible, 
regardless  of  w .  This  is  required  if  we  are  to  make  decision  that  are  applied 
throughout  a  stage.  If  this  restriction  is  inefficient,  shorter  periods  should  be 
used. 

For  example,  suppose  that  Wj  represents  the  current  stage's  unknown 
stream  flow,  and  m*  is  the  current  decision  to  supply  water.  If  we  assume  vv;  is 
zero,  a  manager  is  constrained  to  use  only  previously  stored  water.  Thus, 
regardless  of  the  actual  flow  realized,  the  manager  will  not  make  a  decision  to 
supply  water  that  may  not  be  available. 

If  the  distribution  of  wj  is  unbounded,  we  can  still  apply  this  approach  by 
using  a  sufficiently  unlikely  value  to  bound  the  distribution.  For  example,  in 
flood  control  problems  it  may  be  useful  to  use  the  10,000  year  flood  value.  For  a 
normally  distributed  wj,  a  value  which  has  a  standard  deviation  of  greater  than 
four  often  will  be  sufficiently  conservative  (probability  of  exceedence  <  iff4).  If 
the  control  decisions  are  overly  constrained  by  the  resulting  values  of  wJr  a 
shorter  stage  should  be  used. 
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2.  Types  of  Constraints 

The  type  of  constraints  used  in  a  problem  will  often  be  a  mathematical 
representation  of  physical  limits  of  the  specific  system  being  studied.  Though 
the  constraints  will  depend  on  the  specific  problem,  most  generally  they  will  be 
of  two  forms:  (a)  bounds  on  decision  variables,  u,  and  (b)  bounds  on  state 
variables,  x  and  y.  In  either  case,  the  mathematical  forms  of  the  constraints  are 
reorganized  as  above  to  be  equality  constraints  that  operate  on  u. 

(a)  Bounds  on  Decision  Variables 

These  will  have  the  general  form 

My,  min  —  uj  —  My,max 

where  My  is  single  decision  of  the  entire  set  u.  When  active,  this  constraint  is 
represented  by 

gTu  =  d 

where  g  is  a  vector  of  zeros  except  for  element  k,  which  is  either  one  (if 
d  =  -  Uj  jinn ),  or  minus  one  (if  d  =  My, max)-  Since  the  right-hand  side  is  a  scalar, 
dd/dx  =  0. 

(b)  Bounds  on  State  Variables 

These  will  have  the  general  form 

*y,min  —  Xj  ^  Xy.tnax  and  yy,mm  ^  yj  ^  y/.max 

where  xj  and  yj  are  the  f  th  state  variables  describing  the  current  and  end-of- 
stage  states  of  the  system.  The  constraints  on  these  two  state  variables  will  often 
be  the  same  since  they  represent  the  same  system  characteristic  at  adjacent 
stages.  They  can  be  different,  such  as  when  constraints  change  with  time.  For 
example,  the  maximum  reservoir  storage  available  for  water  supply  may  change 
with  the  season  due  to  capacity  reserved  for  flood  control. 

Constraints  on  x  can  be  ignored  since  we  choose  the  initial  state  of  the 
system,  and  we  only  consider  initial  states  that  are  feasible.  Decisions  affect  y 
through  the  transfer  function  y  =  T(x,u,w)  ;  and,  as  discussed  above,  constraints 
on  y  will  be  reorganized  by  substituting  the  transfer  function  representation 
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and  moving  all  terms  with  u  to  the  left  hand  side.  When  active,  this  constraint 
is  again  represented  by 

gTu  =  d(x) 

and  dd/dx  =  3d/dx . 
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Appendix  B 

Function  Approximation  Through  Hermite  Interpolation 

Given  the  function  value  of  F(x)  and  vector  of  derivatives  dF/dx  at  all 
surrounding  nodes  xt-,  we  can  interpolate  the  values  at  intermediate  states  by 
Hermite  interpolation  following  the  method  presented  by  Kitanidis  [1986]  and 
summarized  by  Foufoula-Georgiou  and  Kitanidis  [1988,  Appendix  A  and  B].  This  is 
accomplished  by  using  a  weighted  sum  of  the  values  at  the  i  =  1,...,2"  corner 
points  of  the  hyper-cube  that  bound  the  intermediate  point  (Figure  B-l).  Unlike 
the  convention  of  the  rest  of  the  text,  here  I  use  subscripts  only  as  identifiers  and 
they  do  not  indicate  derivatives. 

The  weightings  (pi  and  y/ij  applied  to  node  i  of  the  hyper-cube  are 
determined  by  simple  polynomials  that  have  the  following  properties: 


O'th  Order  Value 

1st  Derivatives 

(pi- 

1  at  node  i, 

0  at  all  nodes 

0  otherwise 

w 

0  at  all  nodes  and  all  directions 

1  at  node  i  and  direction  /, 

0  otherwise. 

Given  these  properties  of  (pi  and  y/jj,  the  values  of  the  function  and  its  derivative 
can  be  approximated  as 

2n 

F  =  X  {Ftfi  +  ^fi) 

i-\ 

3F  _  Y  r p dpi  dFj  3y; 

dx  _  “  1  ‘ax  ax  3x  J 

where  \y,-  =  (vfa,...,V'»n)T  • 

The  simplest  polynomials  that  we  can  use  to  specify  (pi  and  y/iy  has  non¬ 
zero  third  derivatives  and  continuous  second  derivatives.  Thus  we  can  define  a 
continuous  second  derivative  of  F(x)  to  be 

a 2f  _  y  fF  a2</>i  BFj  ay, 

3x3x  “  ‘3x3x  dx  dxdx 
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where 


BFj  ay,-  =  y  agj  a2^- 

ax  axax  rf 1  ax  axax * 

;=i 

Let  xa  and  xb  be  vectors  containing  respectively  the  upper  and  lower 

bounds  of  the  hyper-cube  in  each  of  the  k  =  1 dimensions.  The  dimensions 

of  the  hyper-cube  are  then  Ax  =  xa-xb.  The  location  of  the  intermediate  point 

can  then  be  defined  as  %k  =  (jc *  -  x^)IAxk  and  0  <  ^  <  1 . 

The  distance  of  the  point  from  any  given  corner  x,  can  then  be  specified 
T 

as  t\  =  (T}i,...rin)  where 


77*=  &  if  Xk=xkb 

7?*  =(!-&)  if  xk  =  xka  . 

The  weighting  polynomials  can  be  defined  as  follows 

6  =  RP 


Vij  =  77/1  -  TJj)P 


where  7?  =  1  +  77*  -  2^  rjk2 

k=  1  *=1 


p  =  n  n-77*} . 

*=i 


The  above  equations  could  also  use  vector  notation;  thus  R  =  1  +  (1  -  2Ti)Tq. 
The  first  derivatives  of  these  polynomials  are  then 


where 


dQj 

dxs 


dih 

dxs 


[(1  -  4t7j)P  -  RPS] 


dy/ij 

_  977, 

-1 

TP. 

dxs 

a^5 

P  a,b 

n 

■  n 

U- 

■77*) 

*=1 

k*ayb 

T  =  (1  -  77/Xl  -  377/  if  5  =j 
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T  =  Hflty- 1) 


if  s*j 


The  second  derivatives  of  these  polynomials  are 

d2<t>i  _  dry  drjs  c 
-  — 

dxjdxs  dxi  dxs 


where 

$s,s  ~~  6(2  T]s  ■  1  s 

if 

l  =  s 

Su  =  (47],  - 1  )Ps  +  (477/  -  m  +  Mis 

if 

l*s 

and 

oVy  _  377/ 977s 
dxidxs  dxi  dxs  \dxj  J 

where 

Wjjj  =  (6IJ-4 )Pj 

if 

l  —  s  =j 

Wqj  =  OVj-DPi 

if 

l*j,  S=j 

=  (3u,-l)P, 

if 

i  =;'/  s  *i 

=  0 

if 

l*s,  s*j 

Wljj  —  V [/(l  "  Vj)^lyS 

if 

l*s,  l*j,  s*j 

This  presentation  is  essentially  the  same  as  that  of  Foufoula-Georgiou  and 
Kitanidis  [1988]  except  for  correcting  equation  B 7. 
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Hypercube  of  3  Dimensions 


•  discrete  x(i)  =  [xj(i),  x2(i),  x3(i)],  for  i  =  1, 2,  3, ... ,  23 

0  intermediate  x  =  [xj,  x2,  x3] 

or  =  [(x^  +  ^jAxj),  (x2b  +  ^2Ax2),  (x3b  +  ^3Ax3)] 

or  =  [(XjW  +  rijUJAxj),  (x2(1)  +  ti2(1)Ax2),  (x3(1)  +  ti3(1)Ax3)] 

or  =  [(x1(2)  -  ■n1(2)Ax1),  (x2(2)  +  T)2(2)Ax2),  (x3(2)  +  T|3(2)A x3)] 

or  =  [(x1(3)  +  'n1(3)Ax1),  (x2(3)  -  il2(3)Ax2),  (x3(3)  +  T]3(3)Ax3)] 

•  •  • 

or  =  [(xj(8)  -r|1(8)Ax1),  (x2(8)  -  ti2(8)Ax2),  (x3(8)  -  ti3(8)Ax3)] 


Fig.  B-l.  Illustration  of  a  Hypercube  of  3  Dimensions 
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Appendix  C 

Development  of  Derivatives 

The  most  general  form  of  the  objective  function  can  be  represented  as  any 
function  f  =  f{x, u,w,y)  where  u*(x),  w(x) ,  and  y  =  T(x,u,w)  .  Based  on  these 
functional  relationships,  the  following  relationships  can  be  derived: 

dy=5y  du=§H.  dw  _ 

du  9u  '  dx  dx  '  dx  dx  ' 

dy  _  dy  |  dy  du  |  dy  dw 

dx  dx  3u  dx  dw  dx 

In  the  process  of  developing  the  following  derivatives,  I  assume  that 
y  =  T(x,u,w)  is  a  linear  function.  This  allows  us  to  drop  second-order 
derivatives  of  y  and  improve  the  presentation's  clarity.  If  the  transition  function 
were  non-linear,  these  terms  could  be  preserved  with  only  minor  modification  to 
the  following  presentation.  To  the  more  general  formulation,  I  apply  the 
objective  function  proposed  by  Foufoula-Georgiou  and  Kitanidis  [1988]  for  reservoir 
management: 

/=  C(x,u)  +  F(y)  (Cl) 

where  C(x,u)  is  a  single-stage  cost  function,  and  F(y)  is  a  cost  function 
incorporating  expected  costs  of  future  stages. 

A.  Derivatives  of  the  Objective  Function,  /(x,u,w,y) 

(1)  Derivatives  of  objective  function  with  respect  to  decisions: 

V«f  =  /a  +/,  (Cla) 


and,  using  the  more  specific  objective  function  (Cl), 

V/=CU+Fy^ 


(Clb) 


(2)  Second  derivatives  of  objective  function  with  respect  to  decisions: 
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Applying  the  above  relationship  for  V,/  and  dropping  second  order  derivatives 
of  y: 


v  ,v  at  -  f  +  ^ f  +  f  ( 

VutVuO  -  +  8u  /yu  +  /uy  ^  ^  /yy  ^  ,  ' 

t  5y  ^  Sy  / 

and,  more  specifically,  VU(V„/)  =  Cuu  +  ^  ' 

(3)  Derivatives  of  objective  function  with  respect  to  initial  state: 

V-A  +*£• 

Expanding  dy/dx,  we  can  re-express  this  in  terms  of  V,/  as  follows 

and  V/=Vu|Uc,  +  F,(|-  +  J^). 

(4)  Derivatives  of  objective  function  with  respect  to  initial  state  and  decisions 
(Note  that  VU(VX/)T  =  (VX(V1/)T)T  ): 

V..(VAT  =  +  +  aw  +  djT  _ 


(C2a) 


(C2b) 


(C3a) 


(C3b) 


8u  dx  dw  dx 


y  dx 


Applying  the  relationship  for  VJ  and  dropping  second  order  derivatives  of  y: 

VufV^T  =  fm  +/„  57  +  57  (r.  +/»,  §j)  +  57  (a™  +/»j  57) 

dy\ 


.  dyT(f  +f  §y 

dx  r*  /),y  du 


Expanding  dy/dx,  we  can  re-express  this  in  terms  of  Vu(Vu/)T  as  follows 

frK+/*'lr) 

I) 


(C4a) 
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and  Wtf  .  |T Vu(V/7  +  CIU  +  (£  +  £  .  (C4b) 

The  above  form  is  useful  since  terms  with  du/dx  are  aggregated. 

Alternately,  (C3a)  and  (C3b)  can  be  expressed  as 

V  =  C,  +  C„  ^  +  F,g-  (C5a) 

V.(V,fiT  =  C„  +  |H.TC„u  +  ,  (C5b) 

but  these  forms  are  not  used  except  for  comparison  with  the  work  of  Foufoula- 
Georgiou  and  Kitanidis  [1988]. 

MODIFICATION  OF  PRIOR  WORK 

Foufoula-Georgiou  and  Kitanidis  [1988]  proposed  the  following  derivations 


VJ-  =  V„C  +  V,  , 

(C16) 

=  C„U  +  F„^L  , 

(Cl  7) 

V-w.<f*v,.£+£g). 

(C18) 

Vx(V„flT  =  CUX  +  |^T F yy  . 

(C19) 

It  is  useful  to  first  note  the  following  relationships: 

V„c  =  C„  ,  V,F  =  F,  ,  V,(V/)t  =  F„  , 

V„C  =  Cx  +  C„  ,  and  Vv(VuC)T  =  C0<  +  C„„  ^  . 

The  following  observations  are  made  in  comparing  the  results  of  the 
derivations  presented  here  and  those  of  the  original  work: 

(1)  Equation  (C16)  is  the  same  because  of  the  equivalence  of  the  general 
and  partial  derivatives. 
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(2)  Equation  (Cl  7)  is  also  the  same  except  for  the  inconsistent  use  of  the 
symbol  fun  which  implies  a  partial  derivative,  when  the  general  derivative  is 
desired. 

(3)  The  first  two  terms  of  equation  (Cl  8)  are  different  because  of  the 
incorrect  use  of  the  general  derivative  symbols  V  and  d  when  partial 
derivatives  should  have  been  used: 

VXC  =  Cx  +  Cu^-  *  VXC  +  Vu  . 

dx  dx 

(4)  The  definition  for  Vu(Vx/)T  (equation  C19)  in  the  original  work  is 
absent  the  final  term  dy/dx .  Also  the  first  two  terms  are  different  because  of  the 
incorrect  use  of  the  partial  derivative  notation  Cu»  when  the  general  derivative 
should  have  been  specified: 

VX(VUC)T  =  Cux  +  Cuu  0^-  ^  Cux 

(5)  The  total  derivative 

dy  _  dy  dy  du  dy  dw 
dx  dx  +  du  dx  +  dw  dx 

presented  here  is  consistent  with  the  original  work  since,  in  the  absence  of 
autocorrelation,  dw/dx  -  0 ;  thus 

dy  _  dy  |  dy  du 
dx  dx  du  dx 

When  modeling  a  system  with  autocorrelation,  we  need  the  longer  form  to 
evaluate  V/  and  Vu(Vx/)T. 

(6)  Because  of  autocorrelation,  the  matrix  dy/dx  is  not  constant,  and  must 
be  recalculated  for  each  X(,)  and  for  each  W(*) .  Since 

dypt)  _  dy  dy  du  |  dy(^  dw^) 
dx  dx  du  dx  d\V(*)  dx 

each  stochastic  realization  and  discrete  state  affects  the  calculated  value  through 
the  additional  terms  in  this  formulation.  Even  with  a  linear  transition  function, 
the  above  total  derivative  is  not  constant  because  of  the  effect  that  state  variables 
have  on  conditioning  the  distribution  of  w . 
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B.  Approximation  by  Discrete  Realizations 

The  values  derived  from  the  above  formulas  are  uncertain  because  of  the 
presence  of  stochastic  variables  w.  We  can  estimate  these  values  by  using 
probability  weighted  realizations  w(i)  to  approximate  their  expected  values.  We 
can  choose  discrete  realizations  to  span  the  multi-variate  distribution  of  w 
much  as  we  can  chose  states  to  span  x.  The  objective  function  (Cl)  is 

evaluated  as 

/  =  C(x,u)  +  X  {PkF(y(kj))  (C6) 

k 


where  p*  are  the  weights  and  y<*)  =  T(x,u,w^))  .  Similarly, 

v«f  =  c„  +  X  (J't'W  (C7) 

k 

=  C„„  +  |-TX  (C8) 

v/  =  v0y|  +  c,  +  X  WW  g-  ♦  I  Ww  gg  ^f-1 


v»(VxOT  = 


V„(VuflT  +  Cxo+  gTX  g 


v  r  dw(k)Tdy(k)  Tc,  , 

£  to-ar  Fyy’(*)} 


(CIO) 


We  can  calculate  values  for  equations  (C6)  -  (CIO)  quickly  once  we  determine 


k 

(Cll) 

X  {p^y.wl  / 

k 

(C12) 

X  {PkFyyXk) )  / 

k 

(C13) 

V  r  r  ay(/k)  dw(*>  1 

?  (P‘F*<‘>  3w<t,  3x  1  ' 

(C14) 

V  f_  r-  ay(*>  3w(*) , 

f  ,P*  ”■<*>  3w(«  3x  1  ’ 

(C15) 
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C.  Solving  Gradients  du*/dx 

Derivatives  of  the  Lagrangian  formed  from  the  objective  function  and 
constraints  yield 


V„f  +  GTA  =  0 


and  Gu  -  d  =  0 

where  both  u*  and  A  are  implicit  functions  of  x.  Derivatives  with  respect  to  x 
are 


and 

where  the  term  Vu(Vx/)T 
into  two  parts 


V„(Vxf)T  +  GT^  =  0  , 

pdu*  dd  _  o 
dx  dx 

is  defined  in  Appendix  D.  This  term  can  be  divided 


v"<Vs°T  =  +  D 


with  D  -/xu+/x,  +  gj 


ay  awT/  ,  /ay  ay  aw  w 

^-+ —  (/wu+/wyau)  +  [dx  aw  axJru  /] 


yy  au 


This  form  has  the  advantage  that  terms  with  du*/dx  are  grouped  and  permits  us 
to  solve  du*/dx  by  the  simultaneous  equations 


V„(V„0T 

gt 

du*/dx 

-D 

G 

0 

.  dA/dx  . 

Vxd  . 

where  D  is  defined  as  above,  and  Vxd  =  dx  . 

Using  the  objective  function  proposed  by  Foufoula-Georgiou  and  Kitanidis 

[1988], 


D  —  CXu  + 


ay  .  ay  awY T  3y 
t3x  aw  ax  /  yy  au 


Because  D  contains  stochastic  elements,  its  expected  value  is  the  probability 
weighted  estimate 
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D  -  Cxu  +  X  {PkF yy.(*)}  tH,,  +  X 


dy_ 

du 


9w(fc)  T3y(t) 


3y 


dx  dw 


’(*) 


W  Sr 


du 


where  the  terms  involved  in  summation  are  the  same  as  calculated  in  Appendix 
D. 


MODIFIC  ATION  OF  PRIOR  WORK 

The  original  work  of  Foufoula-Georgiou  and  Kitanidis  [1988]  specified  that, 
in  a  deterministic  problem,  the  upper  right-hand-side  sub-matrix  [-  D]  was 

"(Cux  +  !toFyy)  ' 

The  transpose  of  the  matrix  defined  by  this  formula  yields  nearly  the  same  result 
as  the  present  formulation  in  the  absence  of  autocorrelation.  Without  state 
variables  that  incorporate  autocorrelation,  dw/dx  is  a  zero  matrix  and  the  term 
dy/dx  will  often  be  an  identity  matrix. 

Foufoula-Georgiou  and  Kitanidis  [1988]  refer  to  the  above  formula  as 
-  Vx(Vy/)T  which  is  incorrect.  The  correct  value  can  be  found  above. 
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